
      SUBROUTINE DECODE(RHO,C,S)
      IMPLICIT NONE
*     .. Scalar Arguments ..
      DOUBLE PRECISION C,RHO,S
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS,SQRT
*     ..
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      DOUBLE PRECISION TWO
      PARAMETER (TWO=2.0D0)
*     ..
      IF (RHO.EQ.ONE) THEN
          C = ZERO
          S = ONE

      ELSE IF (ABS(RHO).LT.ONE) THEN
          S = TWO*RHO
          C = SQRT(ONE-S**2)

      ELSE
          C = TWO/RHO
          S = SQRT(ONE-C**2)
      END IF

      RETURN

      END
      SUBROUTINE DINIT(N,ALPHA,DX,INCX)
      IMPLICIT NONE

*     Initialises a vector x with a scalar alpha.
*     Modified from dcopy, BLAS Level 1.
*     Rudnei Dias da Cunha, 14/6/93.
*

*     copies a vector, x, to a vector, y.
*     uses unrolled loops for increments equal to one.
*     jack dongarra, linpack, 3/11/78.
*
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA
      INTEGER INCX,N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION DX(*)
*     ..
*     .. Local Scalars ..
      INTEGER I,IX,M,MP1
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC MOD
*     ..
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1) GO TO 20
*
*        code for unequal increments or equal increments
*          not equal to 1
*
      IX = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      DO 10 I = 1,N
          DX(IX) = ALPHA
          IX = IX + INCX
   10 CONTINUE
      RETURN
*
*        code for both increments equal to 1
*
*
*        clean-up loop
*
   20 M = MOD(N,7)
      IF (M.EQ.0) GO TO 40
      DO 30 I = 1,M
          DX(I) = ALPHA
   30 CONTINUE
      IF (N.LT.7) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,7
          DX(I) = ALPHA
          DX(I+1) = ALPHA
          DX(I+2) = ALPHA
          DX(I+3) = ALPHA
          DX(I+4) = ALPHA
          DX(I+5) = ALPHA
          DX(I+6) = ALPHA
   50 CONTINUE
      RETURN

      END

      SUBROUTINE DMACHCONS(WHAT,RESULT)

* These values are for IEEE-754 arithmetic

*     .. Parameters ..
      DOUBLE PRECISION MACHEPS
      PARAMETER (MACHEPS =     1.1102230246252D-16)
      DOUBLE PRECISION UNDERFLOW
      PARAMETER (UNDERFLOW =     2.2250738585072D-307)
      DOUBLE PRECISION OVERFLOW
      PARAMETER (OVERFLOW =     1.7976931348623D+308)
*     ..
*     .. Scalar Arguments ..
      DOUBLE PRECISION RESULT
      CHARACTER WHAT
*     ..
      IF ((WHAT.EQ.'M') .OR. (WHAT.EQ.'m')) THEN
          RESULT = MACHEPS
      ELSE IF ((WHAT.EQ.'U') .OR. (WHAT.EQ.'u')) THEN
          RESULT = UNDERFLOW
      ELSE IF ((WHAT.EQ.'O') .OR. (WHAT.EQ.'o')) THEN
          RESULT = OVERFLOW
      END IF

      RETURN

      END
      DOUBLE PRECISION FUNCTION DSETRHSSTOP(B,WRK,EPSILON,IPAR,PRECONL,
     +                 PDNRM)
      IMPLICIT NONE

*     .. Scalar Arguments ..
      DOUBLE PRECISION EPSILON
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),WRK(*)
      INTEGER IPAR(*)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM
      EXTERNAL PDNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL PRECONL
*     ..
*     .. Local Scalars ..
      INTEGER LOCLEN,STOPTYPE
*     ..
      LOCLEN = IPAR(4)
      STOPTYPE = IPAR(9)
      IF ((STOPTYPE.EQ.1) .OR. (STOPTYPE.EQ.4) .OR.
     +    (STOPTYPE.EQ.7)) THEN
*  ||r||<epsilon or ||Q1r||<epsilon ||x(k)-x(k-1)||<epsilon
          DSETRHSSTOP = EPSILON

      ELSE IF ((STOPTYPE.EQ.2) .OR. (STOPTYPE.EQ.3) .OR.
     +         (STOPTYPE.EQ.5)) THEN
*  ||r||<epsilon||b|| or sqrt(r(Q1r))<epsilon||b|| or
*  ||Q1r||<epsilon||b||
          DSETRHSSTOP = EPSILON*PDNRM(LOCLEN,B,IPAR)

      ELSE IF (STOPTYPE.EQ.6) THEN
*  ||Q1r||<epsilon||Q1b||
          CALL PRECONL(B,WRK,IPAR)
          DSETRHSSTOP = EPSILON*PDNRM(LOCLEN,WRK,IPAR)
      END IF

      RETURN

      END

      SUBROUTINE ENCODE(RHO,C,S)
      IMPLICIT NONE
*     .. Scalar Arguments ..
      DOUBLE PRECISION C,RHO,S
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS,SIGN
*     ..
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      DOUBLE PRECISION TWO
      PARAMETER (TWO=2.0D0)
*     ..
      IF (C.EQ.ZERO) THEN
          RHO = ONE

      ELSE IF (ABS(S).LT.ABS(C)) THEN
          RHO = SIGN(ONE,C)*S/TWO

      ELSE
          RHO = TWO*SIGN(ONE,S)/C
      END IF

      RETURN

      END

      SUBROUTINE GIVENS(A,B,C,S)
      IMPLICIT NONE
*     .. Scalar Arguments ..
      DOUBLE PRECISION A,B,C,S
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION TAU
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS,SQRT
*     ..
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
*     ..
      IF (B.EQ.ZERO) THEN
          C = ONE
          S = ZERO

      ELSE IF (ABS(B).GT.ABS(A)) THEN
          TAU = -A/B
          S = ONE/SQRT(ONE+TAU**2)
          C = S*TAU

      ELSE
          TAU = -B/A
          C = ONE/SQRT(ONE+TAU**2)
          S = C*TAU
      END IF

      RETURN

      END
      SUBROUTINE PIMDBICG(X,B,WRK,IPAR,DPAR,MATVEC,TMATVEC,PRECONL,
     +                    PRECONR,PDSUM,PDNRM,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=2)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM
      EXTERNAL PDNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS,TMATVEC
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION ALPHA,BETA,DELTA,EPSILON,EXITNORM,RHO,RHO0,
     +                 RHSSTOP,XI
      INTEGER BASISDIM,BLKSZ,CNVRTX,IP,IPTILDE,IR,IRTILDE,IS,ITNO,IW,
     +        IXOLD,IZ,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,PROCID,
     +        STATUS,STEPERR,STOPTYPE
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DOTS(2)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT,DSETRHSSTOP
      EXTERNAL DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,PIMDGETPAR,STOPCRIT
*     ..
      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 9999

      END IF

*  Does not need conversion Y=Q2X for residual
      CNVRTX = 0

*  Set indices for mapping local vectors into wrk
      IR = 1
      IRTILDE = IR + LOCLEN
      IPTILDE = IRTILDE + LOCLEN
      IP = IPTILDE + LOCLEN
      IW = IP + LOCLEN
      IZ = IW + LOCLEN
      IS = IZ + LOCLEN
      IXOLD = IS + LOCLEN

*  Set rhs of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PDNRM)

*  1. r=Q1(b-AQ2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL DCOPY(LOCLEN,B,1,WRK(IP),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IP),1)
              CALL PRECONL(WRK(IP),WRK(IR),IPAR)
          END IF

      ELSE
*     r has been set to Qb in the call to dsetrhsstop
          IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)
          END IF

      END IF

*  2. rtilde=ptilde=p=r
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IRTILDE),1)
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IPTILDE),1)
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)

*  3. rho=dot(rtilde,r)
      DOTS(1) = DDOT(LOCLEN,WRK(IRTILDE),1,WRK(IR),1)
      CALL PDSUM(1,DOTS,IPAR)
      RHO = DOTS(1)

*  4. w=Q1AQ2p
      IF (PRECONTYPE.EQ.0) THEN
          CALL MATVEC(WRK(IP),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.1) THEN
          CALL MATVEC(WRK(IP),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.2) THEN
          CALL PRECONR(WRK(IP),WRK(IZ),IPAR)
          CALL MATVEC(WRK(IZ),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.3) THEN
          CALL PRECONR(WRK(IP),WRK(IW),IPAR)
          CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IW),IPAR)
      END IF

*  5. xi=dot(ptilde,w)
      DOTS(1) = DDOT(LOCLEN,WRK(IPTILDE),1,WRK(IW),1)
      CALL PDSUM(1,DOTS,IPAR)
      XI = DOTS(1)

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      DO 20 ITNO = 1,MAXIT

*  6. alpha=rho/xi
          IF (XI.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 6
              GO TO 9999

          END IF

          ALPHA = RHO/XI

*  7. x=x+alpha*p
          CALL DCOPY(LOCLEN,X,1,WRK(IXOLD),1)
          CALL DAXPY(LOCLEN,ALPHA,WRK(IP),1,X,1)

*  8. r=r-alpha*w
          CALL DAXPY(LOCLEN,-ALPHA,WRK(IW),1,WRK(IR),1)

*  9. check stopping criterion
          CALL STOPCRIT(B,WRK(IR),WRK(IZ),X,WRK(IXOLD),WRK(IW),RHSSTOP,
     +                  CNVRTX,EXITNORM,STATUS,IPAR,MATVEC,TMATVEC,
     +                  PRECONR,PDSUM,PDNRM)

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IZ))

          IF (STATUS.EQ.-5) THEN
              STEPERR = 9
              GO TO 9999
          ELSE IF (STATUS.EQ.0) THEN
              GO TO 9999
          END IF
* 10. rtilde=rtilde-alpha*Q1A^{T}Q2ptilde
          IF (PRECONTYPE.EQ.0) THEN
              CALL TMATVEC(WRK(IPTILDE),WRK(IW),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL TMATVEC(WRK(IPTILDE),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IW),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IPTILDE),WRK(IZ),IPAR)
              CALL TMATVEC(WRK(IZ),WRK(IW),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IPTILDE),WRK(IW),IPAR)
              CALL TMATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IW),IPAR)
          END IF

          CALL DAXPY(LOCLEN,-ALPHA,WRK(IW),1,WRK(IRTILDE),1)

* 11. s=Q1AQ2r
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IR),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(WRK(IR),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IR),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IR),WRK(IS),IPAR)
              CALL MATVEC(WRK(IS),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IS),IPAR)
          END IF

* 12. rho=dot(rtilde,r)
          RHO0 = RHO
          DOTS(1) = DDOT(LOCLEN,WRK(IRTILDE),1,WRK(IR),1)

* 13. delta=dot(rtilde,s)
          DOTS(2) = DDOT(LOCLEN,WRK(IRTILDE),1,WRK(IS),1)

*  Accumulate simultaneously partial values
          CALL PDSUM(2,DOTS,IPAR)
          RHO = DOTS(1)
          DELTA = DOTS(2)

* 14. beta=rho/rho0
          IF (RHO0.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 14
              GO TO 9999

          END IF

          BETA = RHO/RHO0

* 15. p=r+beta*p
          CALL DCOPY(LOCLEN,WRK(IP),1,WRK(IZ),1)
          CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IP),1)

* 16. ptilde=rtilde+beta*ptilde
          CALL DCOPY(LOCLEN,WRK(IPTILDE),1,WRK(IZ),1)
          CALL DCOPY(LOCLEN,WRK(IRTILDE),1,WRK(IPTILDE),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IPTILDE),1)

* 17. w=s+beta*w
          CALL DCOPY(LOCLEN,WRK(IW),1,WRK(IZ),1)
          CALL DCOPY(LOCLEN,WRK(IS),1,WRK(IW),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IW),1)

* 18. xi=delta-beta^2*xi
          XI = DELTA - BETA**2*XI

   20 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMDBICGSTAB(X,B,WRK,IPAR,DPAR,MATVEC,PRECONL,PRECONR,
     +                        PDSUM,PDNRM,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=2)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM
      EXTERNAL PDNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION ALPHA,BETA,EPSILON,EXITNORM,KAPPA,MACHEPS,OMEGA,
     +                 RHO,RHO0,RHSSTOP,XI
      INTEGER BASISDIM,BLKSZ,CNVRTX,IP,IQ,IR,IRTILDE,IS,IT,ITNO,IV,IW,
     +        IXOLD,IZ,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,PROCID,
     +        STATUS,STEPERR,STOPTYPE
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DOTS(2)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT,DSETRHSSTOP
      EXTERNAL DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,DINIT,DMACHCONS,PIMDGETPAR,STOPCRIT
*     ..
      CALL DMACHCONS('M',MACHEPS)

      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 9999

      END IF

*  Does not need conversion Y=Q2X for residual
      CNVRTX = 0

*  Set indices for mapping local vectors into wrk
      IR = 1
      IRTILDE = IR + LOCLEN
      IP = IRTILDE + LOCLEN
      IQ = IP + LOCLEN
      IS = IQ + LOCLEN
      IT = IS + LOCLEN
      IV = IT + LOCLEN
      IW = IV + LOCLEN
      IZ = IW + LOCLEN
      IXOLD = IZ + LOCLEN

*  Set rhs of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PDNRM)

*  1. r=Q1(b-AQ2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL DCOPY(LOCLEN,B,1,WRK(IP),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IP),1)
              CALL PRECONL(WRK(IP),WRK(IR),IPAR)
          END IF

      ELSE
*     r has been set to Qb in the call to dsetrhsstop
          IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)
          END IF

      END IF

*  2. rtilde=r
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IRTILDE),1)

*  3. p=v=0
      CALL DINIT(LOCLEN,ZERO,WRK(IP),1)
      CALL DINIT(LOCLEN,ZERO,WRK(IV),1)

*  4. rho=alpha=omega=1
      RHO = ONE
      ALPHA = ONE
      OMEGA = ONE

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      DO 20 ITNO = 1,MAXIT

*  5. rho=dot(rtilde,r)
          RHO0 = RHO
          DOTS(1) = DDOT(LOCLEN,WRK(IRTILDE),1,WRK(IR),1)
          CALL PDSUM(1,DOTS,IPAR)
          RHO = DOTS(1)

*  6. beta=rho*alpha/(rho0*omega)
          KAPPA = RHO0*OMEGA
          IF (KAPPA.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 6
              GO TO 9999

          END IF

          BETA = RHO*ALPHA/KAPPA

*  7. p=r+beta*(p-omega*v)
          CALL DAXPY(LOCLEN,-OMEGA,WRK(IV),1,WRK(IP),1)
          CALL DCOPY(LOCLEN,WRK(IP),1,WRK(IW),1)
          CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IW),1,WRK(IP),1)

*  8. v=Q1AQ2p
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IP),WRK(IV),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(WRK(IP),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IV),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IP),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IV),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IP),WRK(IV),IPAR)
              CALL MATVEC(WRK(IV),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IV),IPAR)
          END IF

*  9. xi=dot(rtilde,v)
          DOTS(1) = DDOT(LOCLEN,WRK(IRTILDE),1,WRK(IV),1)
          CALL PDSUM(1,DOTS,IPAR)
          XI = DOTS(1)

* 10. alpha=rho/xi
          IF (XI.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 10
              GO TO 9999

          END IF

          ALPHA = RHO/XI

* 11. s=r-alpha*v
          CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IS),1)
          CALL DAXPY(LOCLEN,-ALPHA,WRK(IV),1,WRK(IS),1)

* 12. if ||s||<breaktol then soft-breakdown has occurred
          KAPPA = PDNRM(LOCLEN,WRK(IS),IPAR)
          IF (KAPPA.LT.MACHEPS) THEN
              STATUS = -2
              STEPERR = 12
              CALL DAXPY(LOCLEN,ALPHA,WRK(IP),1,X,1)
              GO TO 9999

          END IF

* 13. t=Q1AQ2s
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IS),WRK(IT),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(WRK(IS),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IT),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IS),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IT),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IS),WRK(IT),IPAR)
              CALL MATVEC(WRK(IT),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IT),IPAR)
          END IF

* 14. omega=dot(t,s)/dot(t,t)
          DOTS(1) = DDOT(LOCLEN,WRK(IT),1,WRK(IT),1)
          DOTS(2) = DDOT(LOCLEN,WRK(IT),1,WRK(IS),1)

*  Accumulate simultaneously partial values
          CALL PDSUM(2,DOTS,IPAR)

          IF (DOTS(1).EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 14
              GO TO 9999

          END IF

          OMEGA = DOTS(2)/DOTS(1)

* 15. x=x+alpha*p+omega*s
          CALL DCOPY(LOCLEN,X,1,WRK(IXOLD),1)
          CALL DAXPY(LOCLEN,ALPHA,WRK(IP),1,X,1)
          CALL DAXPY(LOCLEN,OMEGA,WRK(IS),1,X,1)

* 16. r=s-omega*t
          CALL DCOPY(LOCLEN,WRK(IS),1,WRK(IR),1)
          CALL DAXPY(LOCLEN,-OMEGA,WRK(IT),1,WRK(IR),1)

* 17. check stopping criterion
          CALL STOPCRIT(B,WRK(IR),WRK(IZ),X,WRK(IXOLD),WRK(IW),RHSSTOP,
     +                  CNVRTX,EXITNORM,STATUS,IPAR,MATVEC,MATVEC,
     +                  PRECONR,PDSUM,PDNRM)

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IZ))

          IF (STATUS.EQ.-5) THEN
              STEPERR = 17
              GO TO 9999
          ELSE IF (STATUS.EQ.0) THEN
              GO TO 9999
          END IF
   20 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMDCG(X,B,WRK,IPAR,DPAR,MATVEC,PRECONL,PRECONR,PDSUM,
     +                  PDNRM,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=2)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM
      EXTERNAL PDNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION ALPHA,BETA,DELTA,EPSILON,EXITNORM,RDOTR,RDOTR0,
     +                 RHSSTOP,XI
      INTEGER BASISDIM,BLKSZ,CNVRTX,IP,IR,IS,ITNO,IW,IXOLD,IZ,LDA,
     +        LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,PROCID,STATUS,STEPERR,
     +        STOPTYPE
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DOTS(2)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT,DSETRHSSTOP
      EXTERNAL DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,PIMDGETPAR,STOPCRIT
*     ..
      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 9999

      END IF

*  Does not need conversion Y=Q2X for residual
      CNVRTX = 0

*  Set indices for mapping local vectors into wrk
      IR = 1
      IP = IR + LOCLEN
      IW = IP + LOCLEN
      IZ = IW + LOCLEN
      IS = IZ + LOCLEN
      IXOLD = IS + LOCLEN

*  Set rhs of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PDNRM)

*  1. r=Q1(b-AQ2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL DCOPY(LOCLEN,B,1,WRK(IP),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IP),1)
              CALL PRECONL(WRK(IP),WRK(IR),IPAR)
          END IF

      ELSE
*     r has been set to Qb in the call to dsetrhsstop
          IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)
          END IF

      END IF

*  2. p=r
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)

*  3. rdotr=dot(r,r)
      DOTS(1) = DDOT(LOCLEN,WRK(IR),1,WRK(IR),1)
      CALL PDSUM(1,DOTS,IPAR)
      RDOTR = DOTS(1)

*  4. w=Q1AQ2p
      IF (PRECONTYPE.EQ.0) THEN
          CALL MATVEC(WRK(IP),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.1) THEN
          CALL MATVEC(WRK(IP),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.2) THEN
          CALL PRECONR(WRK(IP),WRK(IZ),IPAR)
          CALL MATVEC(WRK(IZ),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.3) THEN
          CALL PRECONR(WRK(IP),WRK(IW),IPAR)
          CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IW),IPAR)
      END IF

*  5. xi=dot(p,w)
      DOTS(1) = DDOT(LOCLEN,WRK(IP),1,WRK(IW),1)
      CALL PDSUM(1,DOTS,IPAR)
      XI = DOTS(1)

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      DO 20 ITNO = 1,MAXIT

*  6. alpha=rdotr/xi
          IF (XI.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 6
              GO TO 9999

          END IF

          ALPHA = RDOTR/XI

*  7. x=x+alpha*p
          CALL DCOPY(LOCLEN,X,1,WRK(IXOLD),1)
          CALL DAXPY(LOCLEN,ALPHA,WRK(IP),1,X,1)

*  8. r=r-alpha*w
          CALL DAXPY(LOCLEN,-ALPHA,WRK(IW),1,WRK(IR),1)

*  9. check stopping criterion
          CALL STOPCRIT(B,WRK(IR),WRK(IZ),X,WRK(IXOLD),WRK(IS),RHSSTOP,
     +                  CNVRTX,EXITNORM,STATUS,IPAR,MATVEC,MATVEC,
     +                  PRECONR,PDSUM,PDNRM)

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IZ))

          IF (STATUS.EQ.-5) THEN
              STEPERR = 9
              GO TO 9999
          ELSE IF (STATUS.EQ.0) THEN
              GO TO 9999
          END IF
* 10. s=Q1AQ2r
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IR),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(WRK(IR),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IR),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IR),WRK(IS),IPAR)
              CALL MATVEC(WRK(IS),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IS),IPAR)
          END IF

* 11. rdotr=dot(r,r)
          RDOTR0 = RDOTR
          DOTS(1) = DDOT(LOCLEN,WRK(IR),1,WRK(IR),1)

* 12. delta=dot(r,s)
          DOTS(2) = DDOT(LOCLEN,WRK(IR),1,WRK(IS),1)

*  Accumulate simultaneously partial values
          CALL PDSUM(2,DOTS,IPAR)
          RDOTR = DOTS(1)
          DELTA = DOTS(2)

* 13. beta=rdotr/rdotr0
          IF (RDOTR0.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 13
              GO TO 9999

          END IF

          BETA = RDOTR/RDOTR0

* 14. p=r+beta*p
          CALL DCOPY(LOCLEN,WRK(IP),1,WRK(IZ),1)
          CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IP),1)

* 15. w=s+beta*w
          CALL DCOPY(LOCLEN,WRK(IW),1,WRK(IZ),1)
          CALL DCOPY(LOCLEN,WRK(IS),1,WRK(IW),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IW),1)

* 16. xi=delta-beta^2*xi
          XI = DELTA - BETA**2*XI

   20 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMDCGEV(X,B,WRK,IPAR,DPAR,MATVEC,PRECONL,PRECONR,
     +                    PDSUM,PDNRM,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      DOUBLE PRECISION TWO
      PARAMETER (TWO=2.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=6)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM
      EXTERNAL PDNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION A1,ALPHA,B1,BETA,C,D,D1,DELTA,E,EPSILON,EXITNORM,
     +                 F,F0,F1,G,H,H0,H1,MACHEPS,MU1,MUN,RDOTR,RDOTR0,
     +                 RHO,RHO0,RHSSTOP,TA,TA1,TB,TB0,XI,XI0
      INTEGER BASISDIM,BLKSZ,CNVRTX,IDT,IET,IP,IR,IRSLT,IS,ITNO,IW,
     +        IXOLD,IZ,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,PROCID,
     +        STATUS,STEPERR,STOPTYPE
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DOTS(2)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION BISECTION,DDOT,DSETRHSSTOP
      EXTERNAL BISECTION,DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,DMACHCONS,PIMDGETPAR,STOPCRIT
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS,MAX,MIN,SQRT
*     ..

*     .. Scalars in Common ..
      DOUBLE PRECISION OLIMIT,OVERFLOW,ULIMIT,UNDERFLOW
*     ..
*     .. Common blocks ..
      COMMON /MACHCONS1/UNDERFLOW,OVERFLOW
      COMMON /MACHCONS2/ULIMIT,OLIMIT
*     ..
      CALL DMACHCONS('M',MACHEPS)
      CALL DMACHCONS('U',UNDERFLOW)
      ULIMIT = SQRT(UNDERFLOW)*TWO
      CALL DMACHCONS('O',OVERFLOW)
      OLIMIT = SQRT(OVERFLOW)/TWO

      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 9999

      END IF

*  Does not need conversion Y=Q2X for residual
      CNVRTX = 0

*  Set indices for mapping local vectors into wrk
      IR = 1
      IP = IR + LOCLEN
      IW = IP + LOCLEN
      IZ = IW + LOCLEN
      IS = IZ + LOCLEN
      IXOLD = IS + LOCLEN
      IDT = IXOLD + LOCLEN
      IET = IDT + MAXIT + 1

*  Set rhs of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PDNRM)

*  1. r=Q1(b-AQ2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL DCOPY(LOCLEN,B,1,WRK(IP),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IP),1)
              CALL PRECONL(WRK(IP),WRK(IR),IPAR)
          END IF

      ELSE
*     r has been set to Qb in the call to dsetrhsstop
          IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)
          END IF

      END IF

*  2. p=r
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)

*  3. rdotr=dot(r,r)
      DOTS(1) = DDOT(LOCLEN,WRK(IR),1,WRK(IR),1)
      CALL PDSUM(1,DOTS,IPAR)
      RDOTR = DOTS(1)

*  4. w=Q1AQ2p
      IF (PRECONTYPE.EQ.0) THEN
          CALL MATVEC(WRK(IP),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.1) THEN
          CALL MATVEC(WRK(IP),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.2) THEN
          CALL PRECONR(WRK(IP),WRK(IZ),IPAR)
          CALL MATVEC(WRK(IZ),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.3) THEN
          CALL PRECONR(WRK(IP),WRK(IW),IPAR)
          CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IW),IPAR)
      END IF

* 5. xi=dot(p,w)
      DOTS(1) = DDOT(LOCLEN,WRK(IP),1,WRK(IW),1)
      CALL PDSUM(1,DOTS,IPAR)
      XI = DOTS(1)

* 6. Initialise data for eigenvalue estimation
      MU1 = ZERO
      MUN = ZERO
      RHO = SQRT(RDOTR)
*  Use rdotr==rho**2
      IF (RDOTR.EQ.ZERO) THEN
          ITNO = 0
          STATUS = -3
          STEPERR = 6
          GO TO 9999

      END IF

      TA1 = XI/RDOTR
      WRK(IDT) = TA1
      TB = ZERO

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      DO 20 ITNO = 1,MAXIT

*  7. alpha=rdotr/xi
          IF (XI.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 7
              GO TO 9999

          END IF

          ALPHA = RDOTR/XI

*  8. x=x+alpha*p
          CALL DCOPY(LOCLEN,X,1,WRK(IXOLD),1)
          CALL DAXPY(LOCLEN,ALPHA,WRK(IP),1,X,1)

*  9. r=r-alpha*w
          CALL DAXPY(LOCLEN,-ALPHA,WRK(IW),1,WRK(IR),1)

* 10. check stopping criterion
          CALL STOPCRIT(B,WRK(IR),WRK(IZ),X,WRK(IXOLD),WRK(IS),RHSSTOP,
     +                  CNVRTX,EXITNORM,STATUS,IPAR,MATVEC,MATVEC,
     +                  PRECONR,PDSUM,PDNRM)

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IZ))

          IF (STATUS.EQ.-5) THEN
              STEPERR = 10
              GO TO 9999
          ELSE IF (STATUS.EQ.0) THEN
              GO TO 9999
          END IF
* 11. s=Q1AQ2r
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IR),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(WRK(IR),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IR),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IR),WRK(IS),IPAR)
              CALL MATVEC(WRK(IS),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IS),IPAR)
          END IF

* 12. rdotr=dot(r,r)
          RDOTR0 = RDOTR
          DOTS(1) = DDOT(LOCLEN,WRK(IR),1,WRK(IR),1)

* 13. delta=dot(r,s)
          DOTS(2) = DDOT(LOCLEN,WRK(IR),1,WRK(IS),1)

*  Accumulate simultaneously partial values
          CALL PDSUM(2,DOTS,IPAR)
          RDOTR = DOTS(1)
          DELTA = DOTS(2)

* 14. beta=rdotr/rdotr0
          IF (RDOTR0.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 14
              GO TO 9999

          END IF

          BETA = RDOTR/RDOTR0

* 15. p=r+beta*p
          CALL DCOPY(LOCLEN,WRK(IP),1,WRK(IZ),1)
          CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IP),1)

* 16. w=s+beta*w
          CALL DCOPY(LOCLEN,WRK(IW),1,WRK(IZ),1)
          CALL DCOPY(LOCLEN,WRK(IS),1,WRK(IW),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IW),1)

* 17. xi=delta-beta^2*xi
          XI0 = XI
          XI = DELTA - BETA**2*XI

* 18. compute estimates of eigenvalues
          TA = TA1
*  Use rdotr==rho**2.0
          RHO0 = RHO
          RHO = SQRT(RDOTR)

          TA1 = (BETA**TWO*XI0+XI)/RDOTR
          TB0 = TB
          TB = -XI0*BETA/ (RHO0*RHO)
          WRK(IDT+ITNO) = TA

*     Store square of subdiagonal elements to save operations in
*     the bisection routine
          WRK(IET+ITNO-1) = TB*TB

          IF (ITNO.EQ.1) THEN
*     Compute roots of the funression for the determinant
*       | a1 b1 |
*       | b1 a2 |
*     i.e., a quadratic equation
*     p(mu)=mu^2-(a1+a2)mu+a1*a2-b1^2 via quadratic formula; take
*     care of possible cancellation of digits
              A1 = TA + TA1
              B1 = TA*TA1 - WRK(IET+ITNO-1)
              IF (A1.LT.ZERO) THEN
                  D1 = SQRT(A1*A1-4.0D0*B1)
                  E = (A1+D1)/TWO
                  F = TWO*B1/ (A1+D1)
                  MU1 = MIN(E,F)
                  MUN = MAX(E,F)

              ELSE IF (A1.GT.ZERO) THEN
                  D1 = SQRT(A1*A1-4.0D0*B1)
                  E = TWO*B1/ (A1-D1)
                  F = (A1-D1)/TWO
                  MU1 = MIN(E,F)
                  MUN = MAX(E,F)

              ELSE
                  D1 = SQRT(-B1)
                  MU1 = -D1
                  MUN = D1
              END IF

*     Set [c,d] to [mu1,mun]
              C = MU1
              D = MUN

          ELSE

*     Update [c,d] via Gerschgorin's theorem
              G = ABS(TB)
              E = ABS(TB0) + G
              F = ABS(TA)
              F0 = F - E
              F1 = F + E
              H = ABS(TA1)
              H0 = H - G
              H1 = H + G
              C = MIN(C,MIN(F0,H0))
              D = MAX(D,MAX(F1,H1))

*     Update [mu1,mun] via bisection
              MU1 = BISECTION(MACHEPS,WRK(IDT),WRK(IET),ITNO+1,C,MU1,
     +              IRSLT)
              IF (IRSLT.EQ.-1) THEN
                  STATUS = -8
                  STEPERR = 18
                  GO TO 9999
              END IF
              IF (IRSLT.EQ.-2) THEN
                  STATUS = -9
                  STEPERR = 18
                  GO TO 9999
              END IF

              MUN = BISECTION(MACHEPS,WRK(IDT),WRK(IET),ITNO+1,MUN,D,
     +              IRSLT)
              IF (IRSLT.EQ.-1) THEN
                  STATUS = -10
                  STEPERR = 18
                  GO TO 9999
              END IF
              IF (IRSLT.EQ.-2) THEN
                  STATUS = -11
                  STEPERR = 18
                  GO TO 9999
              END IF


          END IF

   20 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM
      DPAR(3) = MU1
      DPAR(4) = MUN
      DPAR(5) = 0.0D0
      DPAR(6) = 0.0D0

      RETURN

      END

      DOUBLE PRECISION FUNCTION BISECTION(DMACHEPS,A,B,K,Y,Z,IRSLT)
      IMPLICIT NONE

*  Computes the root of a function (the determinantal function
*  associated with the Lanczos's tridiagonal matrix) using the
*  bisection method.
*  Possible values of IRSLT: 0 - no floating-point errors
*                           -1 - underflow
*                           -2 - overflow
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      DOUBLE PRECISION TWO
      PARAMETER (TWO=2.0D0)
      INTEGER MAXIT
      PARAMETER (MAXIT=50)
*     ..
*     .. Scalar Arguments ..
      DOUBLE PRECISION DMACHEPS,Y,Z
      INTEGER IRSLT,K
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(*),B(*)
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION P2X,P2Y,X
      INTEGER IT
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DETFUN
      EXTERNAL DETFUN
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS
*     ..
*     .. Scalars in Common ..
      DOUBLE PRECISION OVERFLOW,UNDERFLOW
*     ..
*     .. Common blocks ..
      COMMON /MACHCONS1/UNDERFLOW,OVERFLOW
*     ..
      IRSLT = 0
      X = (Y+Z)/TWO
      IF (ABS(Y-Z).LE. (DMACHEPS* (ABS(Y)+ABS(Z)))) THEN
          BISECTION = X
          RETURN
      END IF

* Compute p2(y) once only
      P2Y = DETFUN(A,B,K,Y)
      IF (P2Y.EQ.UNDERFLOW) THEN
          BISECTION = X
          IRSLT = -1
          RETURN
      END IF
      IF (P2Y.EQ.OVERFLOW) THEN
          BISECTION = X
          IRSLT = -2
          RETURN
      END IF

      DO 10 IT = 1,MAXIT
          P2X = DETFUN(A,B,K,X)
          IF (P2X.EQ.UNDERFLOW) THEN
              BISECTION = X
              IRSLT = -1
              RETURN
          END IF
          IF (P2X.EQ.OVERFLOW) THEN
              BISECTION = X
              IRSLT = -2
              RETURN
          END IF

          IF (SIGN(ONE,P2X).NE.SIGN(ONE,P2Y)) THEN
* There is a root between x and y
* p2(y) remains the same for the next iteration
              Z = X

          ELSE IF (ABS(P2X).LT.DMACHEPS) THEN
* Close enought to zero, accept x as root
              BISECTION = X
              RETURN
          ELSE
* There is a root between x and z
* p2(y) becomes p2(x) for next iteration
              Y = X
              P2Y = P2X
          END IF

          IF (ABS(Y-Z).LE. (DMACHEPS* (ABS(Y)+ABS(Z)))) THEN
              BISECTION = X
              RETURN
          END IF
          X = (Y+Z)/TWO
   10 CONTINUE
      BISECTION = X
      RETURN

      END

      DOUBLE PRECISION FUNCTION DETFUN(A,B,K,X)
      IMPLICIT NONE


*  Compute determinant function at x; the squares of
*  the subdiagonal elements of T_{k} are stored in b
*     .. Scalar Arguments ..
      DOUBLE PRECISION X
      INTEGER K
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(*),B(*)
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION P0,P1,P2
      INTEGER I
*     ..
*     .. Scalars in Common ..
      DOUBLE PRECISION OLIMIT,OVERFLOW,ULIMIT,UNDERFLOW
*     ..
*     .. Common blocks ..
      COMMON /MACHCONS1/UNDERFLOW,OVERFLOW
      COMMON /MACHCONS2/ULIMIT,OLIMIT
*     ..
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS
*     ..
      P1 = A(1) - X
      P2 = (A(2)-X)*P1 - B(1)
      DO 10 I = 3,K
          P0 = P1
          P1 = P2
          P2 = (A(I)-X)*P1 - B(I-1)*P0
          IF ((ABS(P2).GT.ZERO) .AND. (ABS(P2).LE.ULIMIT)) THEN
              DETFUN = UNDERFLOW
              RETURN
          END IF
          IF (ABS(P2).GE.OLIMIT) THEN
              DETFUN = OVERFLOW
              RETURN
          END IF

   10 CONTINUE

      DETFUN = P2
      RETURN

      END
      SUBROUTINE PIMDCGNE(X,B,WRK,IPAR,DPAR,MATVEC,TMATVEC,PRECONL,
     +                    PRECONR,PDSUM,PDNRM,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=2)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM
      EXTERNAL PDNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS,TMATVEC
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION ALPHA,BETA,DELTA,EPSILON,EXITNORM,RDOTR,RDOTR0,
     +                 RHSSTOP,XI
      INTEGER BASISDIM,BLKSZ,CNVRTX,IP,IR,IS,ITNO,IW,IXOLD,IZ,LDA,
     +        LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,PROCID,STATUS,STEPERR,
     +        STOPTYPE
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DOTS(2)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT,DSETRHSSTOP
      EXTERNAL DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,PIMDGETPAR,STOPCRIT
*     ..
      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 9999

      END IF

*  Needs conversion Y=Q2X for residual
      CNVRTX = 1

*  Set indices for mapping local vectors into wrk
      IR = 1
      IP = IR + LOCLEN
      IW = IP + LOCLEN
      IZ = IW + LOCLEN
      IS = IZ + LOCLEN
      IXOLD = IS + LOCLEN

*  Set rhs of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PDNRM)

*  1. r=Q1(b-AA^{T}Q2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL TMATVEC(X,WRK(IP),IPAR)
              CALL MATVEC(WRK(IP),WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-AA^{T})
              CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL TMATVEC(X,WRK(IP),IPAR)
              CALL MATVEC(WRK(IP),WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AA^{T}Q2x
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL PRECONR(X,WRK(IP),IPAR)
              CALL TMATVEC(WRK(IP),WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AA^{T}Q2x)
              CALL DCOPY(LOCLEN,B,1,WRK(IP),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL TMATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IP),1)
              CALL PRECONL(WRK(IP),WRK(IR),IPAR)
          END IF

      ELSE
*     r has been set to Qb in the call to dsetrhsstop
          IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-AA^{T}x)
              CALL TMATVEC(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IP),IPAR)
              CALL PRECONL(WRK(IP),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AA^{T}Q2x)
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL TMATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IR),1)
          END IF

      END IF

*  2. p=r
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)

*  3. rdotr=dot(r,r)
      DOTS(1) = DDOT(LOCLEN,WRK(IR),1,WRK(IR),1)
      CALL PDSUM(1,DOTS,IPAR)
      RDOTR = DOTS(1)

*  4. w=Q1AA^{T}Q2p
      IF (PRECONTYPE.EQ.0) THEN
          CALL TMATVEC(WRK(IP),WRK(IZ),IPAR)
          CALL MATVEC(WRK(IZ),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.1) THEN
          CALL TMATVEC(WRK(IP),WRK(IW),IPAR)
          CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.2) THEN
          CALL PRECONR(WRK(IP),WRK(IW),IPAR)
          CALL TMATVEC(WRK(IW),WRK(IZ),IPAR)
          CALL MATVEC(WRK(IZ),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.3) THEN
          CALL PRECONR(WRK(IP),WRK(IZ),IPAR)
          CALL TMATVEC(WRK(IZ),WRK(IW),IPAR)
          CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IW),IPAR)
      END IF

*  5. xi=dot(p,w)
      DOTS(1) = DDOT(LOCLEN,WRK(IP),1,WRK(IW),1)
      CALL PDSUM(1,DOTS,IPAR)
      XI = DOTS(1)

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      DO 20 ITNO = 1,MAXIT

*  6. alpha=rdotr/xi
          IF (XI.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 6
              GO TO 9999

          END IF

          ALPHA = RDOTR/XI

*  7. x=x+alpha*p
          CALL DCOPY(LOCLEN,X,1,WRK(IXOLD),1)
          CALL DAXPY(LOCLEN,ALPHA,WRK(IP),1,X,1)

*  8. r=r-alpha*w
          CALL DAXPY(LOCLEN,-ALPHA,WRK(IW),1,WRK(IR),1)

*  9. check stopping criterion
          CALL STOPCRIT(B,WRK(IR),WRK(IZ),X,WRK(IXOLD),WRK(IS),RHSSTOP,
     +                  CNVRTX,EXITNORM,STATUS,IPAR,MATVEC,TMATVEC,
     +                  PRECONR,PDSUM,PDNRM)

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IZ))

          IF (STATUS.EQ.-5) THEN
              STEPERR = 9
              GO TO 9999
          ELSE IF (STATUS.EQ.0) THEN
              GO TO 9999
          END IF
* 10. s=Q1AA^{T}Q2r
          IF (PRECONTYPE.EQ.0) THEN
              CALL TMATVEC(WRK(IR),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL TMATVEC(WRK(IR),WRK(IS),IPAR)
              CALL MATVEC(WRK(IS),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IR),WRK(IS),IPAR)
              CALL TMATVEC(WRK(IS),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IR),WRK(IZ),IPAR)
              CALL TMATVEC(WRK(IZ),WRK(IS),IPAR)
              CALL MATVEC(WRK(IS),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IS),IPAR)
          END IF

* 11. rdotr=dot(r,r)
          RDOTR0 = RDOTR
          DOTS(1) = DDOT(LOCLEN,WRK(IR),1,WRK(IR),1)

* 12. delta=dot(r,s)
          DOTS(2) = DDOT(LOCLEN,WRK(IR),1,WRK(IS),1)

*  Accumulate simultaneously partial values
          CALL PDSUM(2,DOTS,IPAR)
          RDOTR = DOTS(1)
          DELTA = DOTS(2)

* 13. beta=rdotr/rdotr0
          IF (RDOTR0.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 10
              GO TO 9999

          END IF

          BETA = RDOTR/RDOTR0

* 14. p=r+beta*p
          CALL DCOPY(LOCLEN,WRK(IP),1,WRK(IZ),1)
          CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IP),1)

* 15. w=s+beta*w
          CALL DCOPY(LOCLEN,WRK(IW),1,WRK(IZ),1)
          CALL DCOPY(LOCLEN,WRK(IS),1,WRK(IW),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IW),1)

* 16. xi=delta-beta^2*xi
          XI = DELTA - BETA**2*XI

   20 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),WRK(IXOLD),IPAR)
          CALL TMATVEC(WRK(IXOLD),X,IPAR)
      ELSE
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL TMATVEC(WRK(IZ),X,IPAR)
      END IF


*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMDCGNR(X,B,WRK,IPAR,DPAR,MATVEC,TMATVEC,PRECONL,
     +                    PRECONR,PDSUM,PDNRM,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=2)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM
      EXTERNAL PDNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS,TMATVEC
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION ALPHA,BETA,DELTA,EPSILON,EXITNORM,RDOTR,RDOTR0,
     +                 RHSSTOP,XI
      INTEGER BASISDIM,BLKSZ,CNVRTX,IP,IR,IS,ITNO,IW,IXOLD,IZ,LDA,
     +        LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,PROCID,STATUS,STEPERR,
     +        STOPTYPE
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DOTS(2)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT,DSETRHSSTOP
      EXTERNAL DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,PIMDGETPAR,STOPCRIT
*     ..
      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 9999

      END IF

*  Does not need conversion Y=Q2X for residual
      CNVRTX = 0

*  Set indices for mapping local vectors into wrk
      IR = 1
      IP = IR + LOCLEN
      IW = IP + LOCLEN
      IZ = IW + LOCLEN
      IS = IZ + LOCLEN
      IXOLD = IS + LOCLEN

*  Set rhs of stopping criteria
*  On (P)CGNR, the rhs vector is actually A^{T}b
      CALL TMATVEC(B,WRK(IZ),IPAR)
      RHSSTOP = DSETRHSSTOP(WRK(IZ),WRK(IR),EPSILON,IPAR,PRECONL,PDNRM)

*  1. r=Q1(A^{T}b-A^{T}AQ2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL MATVEC(X,WRK(IP),IPAR)
              CALL TMATVEC(WRK(IP),WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL DCOPY(LOCLEN,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(A^{T}b-A^{T}Ax)
              CALL MATVEC(X,WRK(IP),IPAR)
              CALL TMATVEC(WRK(IP),WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=A^{T}b-A^{T}AQ2x
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IP),IPAR)
              CALL TMATVEC(WRK(IP),WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL DCOPY(LOCLEN,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(A^{T}b-A^{T}AQ2x)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IP),IPAR)
              CALL TMATVEC(WRK(IP),WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IR),IPAR)
          END IF

      ELSE
*     r has been set to Q(A^{T}b) in the call to dsetrhsstop
          IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(A^{T}b-A^{T}Ax)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL TMATVEC(WRK(IW),WRK(IP),IPAR)
              CALL PRECONL(WRK(IP),WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(A^{T}b-A^{T}AQ2x)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IP),IPAR)
              CALL TMATVEC(WRK(IP),WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IP),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IP),1,WRK(IR),1)
          END IF

      END IF

*  2. p=r
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)

*  3. rdotr=dot(r,r)
      DOTS(1) = DDOT(LOCLEN,WRK(IR),1,WRK(IR),1)
      CALL PDSUM(1,DOTS,IPAR)
      RDOTR = DOTS(1)

*  4. w=Q1A^{T}AQ2p
      IF (PRECONTYPE.EQ.0) THEN
          CALL MATVEC(WRK(IP),WRK(IZ),IPAR)
          CALL TMATVEC(WRK(IZ),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.1) THEN
          CALL MATVEC(WRK(IP),WRK(IW),IPAR)
          CALL TMATVEC(WRK(IW),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.2) THEN
          CALL PRECONR(WRK(IP),WRK(IW),IPAR)
          CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
          CALL TMATVEC(WRK(IZ),WRK(IW),IPAR)

      ELSE IF (PRECONTYPE.EQ.3) THEN
          CALL PRECONR(WRK(IP),WRK(IZ),IPAR)
          CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
          CALL TMATVEC(WRK(IW),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IW),IPAR)
      END IF

*  5. xi=dot(p,w)
      DOTS(1) = DDOT(LOCLEN,WRK(IP),1,WRK(IW),1)
      CALL PDSUM(1,DOTS,IPAR)
      XI = DOTS(1)

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      DO 20 ITNO = 1,MAXIT

*  6. alpha=rdotr/xi
          IF (XI.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 6
              GO TO 9999

          END IF

          ALPHA = RDOTR/XI

*  7. x=x+alpha*p
          CALL DCOPY(LOCLEN,X,1,WRK(IXOLD),1)
          CALL DAXPY(LOCLEN,ALPHA,WRK(IP),1,X,1)

*  8. r=r-alpha*w
          CALL DAXPY(LOCLEN,-ALPHA,WRK(IW),1,WRK(IR),1)

*  9. check stopping criterion
          CALL STOPCRIT(B,WRK(IR),WRK(IZ),X,WRK(IXOLD),WRK(IS),RHSSTOP,
     +                  CNVRTX,EXITNORM,STATUS,IPAR,MATVEC,TMATVEC,
     +                  PRECONR,PDSUM,PDNRM)

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IZ))

          IF (STATUS.EQ.-5) THEN
              STEPERR = 9
              GO TO 9999
          ELSE IF (STATUS.EQ.0) THEN
              GO TO 9999
          END IF
* 10. s=Q1A^{T}AQ2r
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IR),WRK(IZ),IPAR)
              CALL TMATVEC(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(WRK(IR),WRK(IS),IPAR)
              CALL TMATVEC(WRK(IS),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IR),WRK(IS),IPAR)
              CALL MATVEC(WRK(IS),WRK(IZ),IPAR)
              CALL TMATVEC(WRK(IZ),WRK(IS),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IR),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IS),IPAR)
              CALL TMATVEC(WRK(IS),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IS),IPAR)
          END IF

* 11. rdotr=dot(r,r)
          RDOTR0 = RDOTR
          DOTS(1) = DDOT(LOCLEN,WRK(IR),1,WRK(IR),1)

* 12. delta=dot(r,s)
          DOTS(2) = DDOT(LOCLEN,WRK(IR),1,WRK(IS),1)

*  Accumulate simultaneously partial values
          CALL PDSUM(2,DOTS,IPAR)
          RDOTR = DOTS(1)
          DELTA = DOTS(2)

* 13. beta=rdotr/rdotr0
          IF (RDOTR0.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 13
              GO TO 9999

          END IF

          BETA = RDOTR/RDOTR0

* 14. p=r+beta*p
          CALL DCOPY(LOCLEN,WRK(IP),1,WRK(IZ),1)
          CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IP),1)

* 15. w=s+beta*w
          CALL DCOPY(LOCLEN,WRK(IW),1,WRK(IZ),1)
          CALL DCOPY(LOCLEN,WRK(IS),1,WRK(IW),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IZ),1,WRK(IW),1)

* 16. xi=delta-beta^2*xi
          XI = DELTA - BETA**2*XI

   20 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMDCGS(X,B,WRK,IPAR,DPAR,MATVEC,PRECONL,PRECONR,PDSUM,
     +                   PDNRM,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=2)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM
      EXTERNAL PDNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION ALPHA,BETA,EPSILON,EXITNORM,RHO,RHO0,RHSSTOP,XI
      INTEGER BASISDIM,BLKSZ,CNVRTX,IP,IQ,IR,IRTILDE,IS,IT,ITNO,IU,IW,
     +        IXOLD,IZ,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,PROCID,
     +        STATUS,STEPERR,STOPTYPE
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DOTS(1)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT,DSETRHSSTOP
      EXTERNAL DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,PIMDGETPAR,STOPCRIT
*     ..
      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 9999

      END IF

*  Does not need conversion Y=Q2X for residual
      CNVRTX = 0

*  Set indices for mapping local vectors into wrk
      IR = 1
      IRTILDE = IR + LOCLEN
      IP = IRTILDE + LOCLEN
      IQ = IP + LOCLEN
      IS = IQ + LOCLEN
      IT = IS + LOCLEN
      IU = IT + LOCLEN
      IW = IU + LOCLEN
      IZ = IW + LOCLEN
      IXOLD = IZ + LOCLEN

*  Set rhs of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PDNRM)

*  1. r=Q1(b-AQ2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL DCOPY(LOCLEN,B,1,WRK(IP),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IP),1)
              CALL PRECONL(WRK(IP),WRK(IR),IPAR)
          END IF

      ELSE
*     r has been set to Qb in the call to dsetrhsstop
          IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)
          END IF

      END IF

*  2. p=s=rtilde=r
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IRTILDE),1)
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IS),1)

*  3. rho=dot(rtilde,r)
      DOTS(1) = DDOT(LOCLEN,WRK(IRTILDE),1,WRK(IR),1)
      CALL PDSUM(1,DOTS,IPAR)
      RHO = DOTS(1)

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      DO 20 ITNO = 1,MAXIT

*  4. w=Q1AQ2p
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IP),WRK(IW),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(WRK(IP),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IW),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IP),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IP),WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IW),IPAR)
          END IF

*  5. xi=dot(rtilde,w)
          DOTS(1) = DDOT(LOCLEN,WRK(IRTILDE),1,WRK(IW),1)
          CALL PDSUM(1,DOTS,IPAR)
          XI = DOTS(1)

*  6. alpha=rho/xi
          IF (XI.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 6
              GO TO 9999

          END IF

          ALPHA = RHO/XI

*  7. t=s-alpha*w
          CALL DCOPY(LOCLEN,WRK(IS),1,WRK(IT),1)
          CALL DAXPY(LOCLEN,-ALPHA,WRK(IW),1,WRK(IT),1)

*  8. w=s+t
          CALL DCOPY(LOCLEN,WRK(IS),1,WRK(IW),1)
          CALL DAXPY(LOCLEN,ONE,WRK(IT),1,WRK(IW),1)

*  9. x=x+alpha*w
          CALL DCOPY(LOCLEN,X,1,WRK(IXOLD),1)
          CALL DAXPY(LOCLEN,ALPHA,WRK(IW),1,X,1)

* 10. r=r-alpha*Q1AQ2w
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IW),WRK(IU),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IU),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IW),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IU),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IW),WRK(IU),IPAR)
              CALL MATVEC(WRK(IU),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IU),IPAR)
          END IF

          CALL DAXPY(LOCLEN,-ALPHA,WRK(IU),1,WRK(IR),1)

* 11. check stopping criterion
          CALL STOPCRIT(B,WRK(IR),WRK(IZ),X,WRK(IXOLD),WRK(IU),RHSSTOP,
     +                  CNVRTX,EXITNORM,STATUS,IPAR,MATVEC,MATVEC,
     +                  PRECONR,PDSUM,PDNRM)

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IZ))

          IF (STATUS.EQ.-5) THEN
              STEPERR = 11
              GO TO 9999
          ELSE IF (STATUS.EQ.0) THEN
              GO TO 9999
          END IF
* 12. rho=dot(rtilde0,r)
          RHO0 = RHO
          DOTS(1) = DDOT(LOCLEN,WRK(IRTILDE),1,WRK(IR),1)
          CALL PDSUM(1,DOTS,IPAR)
          RHO = DOTS(1)

* 13. beta=rho/rho0
          IF (RHO0.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 13
              GO TO 9999

          END IF

          BETA = RHO/RHO0

* 14. s=r+beta*t
          CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IS),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IT),1,WRK(IS),1)

* 15. w=t+beta*p
          CALL DCOPY(LOCLEN,WRK(IT),1,WRK(IW),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IP),1,WRK(IW),1)

* 16. p=s+beta*w
          CALL DCOPY(LOCLEN,WRK(IS),1,WRK(IP),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IW),1,WRK(IP),1)

   20 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMDCHEBYSHEV(X,B,WRK,IPAR,DPAR,MATVEC,PRECONL,PRECONR,
     +                         PDSUM,PDNRM,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      DOUBLE PRECISION TWO
      PARAMETER (TWO=2.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=6)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM
      EXTERNAL PDNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION AXISISQ,AXISRSQ,D,DELTA,EPSILON,EXITNORM,GAMMA,
     +                 LENGTHI,LENGTHR,RHO,RHSSTOP,SIGMA,SIGMASQ
      INTEGER BASISDIM,BLKSZ,CNVRTX,IK,IR,ITNO,IW,IXOLD,IZ,LDA,LOCLEN,
     +        MAXIT,N,NPROCS,PRECONTYPE,PROCID,STATUS,STEPERR,STOPTYPE
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DSETRHSSTOP
      EXTERNAL DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,DSCAL,DSWAP,PIMDGETPAR,STOPCRIT
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC MAX
*     ..
      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of stop types
      IF ((STOPTYPE.NE.1) .AND. (STOPTYPE.NE.2) .AND.
     +    (STOPTYPE.NE.7)) THEN
          ITNO = 0
          STATUS = -6
          STEPERR = 0
          GO TO 9999

      END IF

*  Does not need conversion Y=Q2X for residual
      CNVRTX = 0

*  Set indices for mapping local vectors into wrk
      IW = 1
      IK = IW + LOCLEN
      IZ = IK + LOCLEN
      IR = IZ + LOCLEN
      IXOLD = IR + LOCLEN

*  Set rhs of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PDNRM)

*  1. Set parameters for iteration
      IF ((DPAR(3).EQ.ZERO) .AND. (DPAR(4).EQ.ZERO) .AND.
     +    (DPAR(5).EQ.ZERO) .AND. (DPAR(6).EQ.ZERO)) THEN
          STATUS = -7
          STEPERR = 1
          GO TO 9999

      ELSE IF (DPAR(5).EQ.DPAR(6)) THEN
*     Eigenvalues are contained in the interval [DPAR(3),DPAR(4)] on
*     the real axis:
*         sigma=(dpar(4)-dpar(3))/(2-dpar(4)-dpar(3))
*         gamma=2/(2-dpar(4)-dpar(3))
          SIGMA = (DPAR(4)-DPAR(3))/ (TWO-DPAR(4)-DPAR(3))
          SIGMASQ = SIGMA*SIGMA
          GAMMA = TWO/ (TWO-DPAR(4)-DPAR(3))

      ELSE IF (DPAR(3).EQ.DPAR(4)) THEN
*     Eigenvalues are contained in the interval [DPAR(5),DPAR(6)] on
*     the imaginary axis:
*         sigma^2=-max(dpar(5),dpar(6))
*         gamma=1
          SIGMASQ = -MAX(DPAR(5),DPAR(6))
          GAMMA = ONE

      ELSE
*     Eigenvalues are complex and contained in the box
*     DPAR(3)<= Real(e) <= DPAR(4) and DPAR(5)<= Imag(e) <= DPAR(6).
*     Compute the minimum bounding ellipse that circumscribes the box;
*     this is defined by its axes a=sqrt(2)*(dpar(4)-dpar(3))/2 (along
*     the real axis) and b=sqrt(2)*(dpar(6)-dpar(5))/2 (along the
*     imaginary axis). The center of the ellipse is d.
*         sigma^2=(a^2+b^2)/(1-d)^2
*         gamma=1/(1-d)
          LENGTHR = (DPAR(4)-DPAR(3))/TWO
          LENGTHI = (DPAR(6)-DPAR(5))/TWO
          AXISRSQ = LENGTHR*LENGTHR*TWO
          AXISISQ = LENGTHI*LENGTHI*TWO
          D = (DPAR(6)+DPAR(5))/TWO
          SIGMASQ = (AXISRSQ-AXISISQ)/ (ONE-D)**2
          GAMMA = ONE/ (ONE-D)

      END IF

*  2. k=gamma*Q1b
      IF (PRECONTYPE.EQ.0) THEN
          CALL DCOPY(LOCLEN,B,1,WRK(IK),1)
          CALL DSCAL(LOCLEN,GAMMA,WRK(IK),1)

      ELSE IF ((PRECONTYPE.EQ.1) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL PRECONL(B,WRK(IK),IPAR)
          CALL DSCAL(LOCLEN,GAMMA,WRK(IK),1)
      END IF

*    xold=x
      CALL DCOPY(LOCLEN,X,1,WRK(IXOLD),1)

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      DO 10 ITNO = 1,MAXIT

*  3. rho
          IF (ITNO.EQ.1) THEN
              RHO = ONE

          ELSE IF (ITNO.EQ.2) THEN
              RHO = ONE/ (ONE-SIGMASQ/TWO)

          ELSE
              RHO = ONE/ (ONE-RHO*SIGMASQ/4.0D0)
          END IF

*  4. w=(I-Q1AQ2)x
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(X,WRK(IZ),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
          END IF

          CALL DCOPY(LOCLEN,X,1,WRK(IW),1)
          CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IW),1)

*  5. x=rho*(gamma*((I-Q1A)x+Q1b)+(1-gamma)*x)+(1-rho)*xold
          DELTA = RHO*GAMMA
          CALL DSCAL(LOCLEN,ONE-RHO,WRK(IXOLD),1)
          CALL DAXPY(LOCLEN,RHO,WRK(IK),1,WRK(IXOLD),1)
          CALL DAXPY(LOCLEN,RHO-DELTA,X,1,WRK(IXOLD),1)
          CALL DAXPY(LOCLEN,DELTA,WRK(IW),1,WRK(IXOLD),1)
          CALL DSWAP(LOCLEN,WRK(IXOLD),1,X,1)

*  6. check stopping criterion
          CALL STOPCRIT(B,WRK(IZ),WRK(IR),X,WRK(IXOLD),WRK(IW),RHSSTOP,
     +                  CNVRTX,EXITNORM,STATUS,IPAR,MATVEC,MATVEC,
     +                  PRECONR,PDSUM,PDNRM)

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IR))

          IF (STATUS.EQ.-5) THEN
              STEPERR = 6
              GO TO 9999
          ELSE IF (STATUS.EQ.0) THEN
              GO TO 9999
          END IF

   10 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,
     +                      NPROCS,PROCID,PRECONTYPE,STOPTYPE,MAXIT,
     +                      ITNO,STATUS,STEPERR,EPSILON,EXITNORM)

c           PIM -- The Parallel Iterative Methods package
c           ---------------------------------------------
c
c                      Rudnei Dias da Cunha
c               Centro de Processamento de Dados,
c         Universidade Federal do Rio Grande do Sul, Brasil
c                              and
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c                          Tim Hopkins
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c ----------------------------------------------------------------------
c
c  Description of parameter arrays
c   IPAR (INPUT)  : integer
c     ipar( 1): lda    (Leading dimension of a)
c           2 : n      (Number of rows/columns of a)
c           3 : blksz  (Size of block of data; used when data is
c                       partitioned using cyclic mode)
c           4 : loclen (Number of elements stored locally;
c                       *PARALLEL: Equal to at least m/nprocs or
c                                  n/procs depending if row or
c                                  column partitioning is used or,
c                                  in the case of cyclic partitioning,
c                                  it is a multiple of either
c                                  m/(blksz*nprocs) or n/(blksz*nprocs).
c                       *SEQUENTIAL: equal to n)
c           5 : basisdim (Dimension of orthogonal basis, used in
c                       GMRES)
c           6 : nprocs (Number of processors)
c           7 : procid (Processor identification)
c           8 : precontype (Type of preconditioning; one of
c                           0 : No preconditioning,
c                           1 : Left preconditioning,
c                           2 : Right preconditioning,
c                           3 : Symmetric preconditioning)
c           9 : stoptype (Type of stopping criteria used)
c          10 : maxit  (Maximum number of iterations allowed)
c
c   IPAR (OUTPUT) : integer
c     ipar(11): itno   (Number of iterations executed)
c          12 : status (On exit of iterative method, one of
c                        0: converged to solution
c                       -1: no convergence has been achieved
c                       -2: "soft"-breakdown, solution may have
c                           been found
c                       -3: "hard"-breakdown, no solution)
c                       -4: conflict in preconditioner and stopping
c                           criterion selected
c                       -5: error in stopping criterion 3, r^{T}z<0)
c          13 : steperr (If status is either -2 or -3, it gives
c                         the step number of the respective algorithm
c                         where a breakdown has occurred. Refer to the
c                         User's Guide for further information)
c
c   RPAR/DPAR (INPUT)  : real/double precision
c     dpar( 1): epsilon (The value of epsilon for use in the
c                        stopping criterion)
c
c   RPAR/DPAR (OUTPUT) : real/double precision
c     dpar( 2): exitnorm (The value of a norm of either the residual
c                         vector or the difference between two 
c                         successive solution estimates according to
c                         the value of stoptype)
c           3,
c           4 : smallest and largest eigenvalues of Q1AQ2 (in the
c               symmetric case) OR smallest and largest real parts
c               (in the nonsymmetric case)
c           5,
c           6 : smallest and largest imaginary parts (only in the
c               nonsymmetric case)
c

C     .. Parameters ..
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=6)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION EPSILON,EXITNORM
      INTEGER BASISDIM,BLKSZ,ITNO,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,
     +        PROCID,STATUS,STEPERR,STOPTYPE
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION DPAR(DPARSIZ)
      INTEGER IPAR(IPARSIZ)
C     ..
      LDA = IPAR(1)
      N = IPAR(2)
      BLKSZ = IPAR(3)
      LOCLEN = IPAR(4)
      BASISDIM = IPAR(5)
      NPROCS = IPAR(6)
      PROCID = IPAR(7)
      PRECONTYPE = IPAR(8)
      STOPTYPE = IPAR(9)
      MAXIT = IPAR(10)
      ITNO = IPAR(11)
      STATUS = IPAR(12)
      STEPERR = IPAR(13)

      EPSILON = DPAR(1)
      EXITNORM = DPAR(2)

      RETURN

      END
      SUBROUTINE PIMDPRTPAR(IPAR,DPAR)

c           PIM -- The Parallel Iterative Methods package
c           ---------------------------------------------
c
c                      Rudnei Dias da Cunha
c               Centro de Processamento de Dados,
c         Universidade Federal do Rio Grande do Sul, Brasil
c                              and
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c                          Tim Hopkins
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c ----------------------------------------------------------------------
c
c  Description of parameter arrays
c   IPAR (INPUT)  : integer
c     ipar( 1): lda    (Leading dimension of a)
c           2 : n      (Number of rows/columns of a)
c           3 : blksz  (Size of block of data; used when data is
c                       partitioned using cyclic mode)
c           4 : loclen (Number of elements stored locally;
c                       *PARALLEL: Equal to at least m/nprocs or
c                                  n/procs depending if row or
c                                  column partitioning is used or,
c                                  in the case of cyclic partitioning,
c                                  it is a multiple of either
c                                  m/(blksz*nprocs) or n/(blksz*nprocs).
c                       *SEQUENTIAL: equal to n)
c           5 : basisdim (Dimension of orthogonal basis, used in
c                       GMRES)
c           6 : nprocs (Number of processors)
c           7 : procid (Processor identification)
c           8 : precontype (Type of preconditioning; one of
c                           0 : No preconditioning,
c                           1 : Left preconditioning,
c                           2 : Right preconditioning,
c                           3 : Symmetric preconditioning)
c           9 : stoptype (Type of stopping criteria used)
c          10 : maxit  (Maximum number of iterations allowed)
c
c   IPAR (OUTPUT) : integer
c     ipar(11): itno   (Number of iterations executed)
c          12 : status (On exit of iterative method, one of
c                        0: converged to solution
c                       -1: no convergence has been achieved
c                       -2: "soft"-breakdown, solution may have
c                           been found
c                       -3: "hard"-breakdown, no solution)
c                       -4: conflict in preconditioner and stopping
c                           criterion selected
c                       -5: error in stopping criterion 3, r^{T}z<0)
c          13 : steperr (If status is either -2 or -3, it gives
c                         the step number of the respective algorithm
c                         where a breakdown has occurred. Refer to the
c                         User's Guide for further information)
c
c   RPAR/DPAR (INPUT)  : real/double precision
c     dpar( 1): epsilon (The value of epsilon for use in the
c                        stopping criterion)
c
c   RPAR/DPAR (OUTPUT) : real/double precision
c     dpar( 2): exitnorm (The value of a norm of either the residual
c                         vector or the difference between two 
c                         successive solution estimates according to 
c                         the value of stoptype)
c           3,
c           4 : smallest and largest eigenvalues of Q1AQ2 (in the
c               symmetric case) OR smallest and largest real parts
c               (in the nonsymmetric case)
c           5,
c           6 : smallest and largest imaginary parts (only in the
c               nonsymmetric case)
c

C     .. Parameters ..
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=6)
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION DPAR(DPARSIZ)
      INTEGER IPAR(IPARSIZ)
C     ..
C     .. Local Scalars ..
      INTEGER I
C     ..
      WRITE (6,FMT=10) (IPAR(I),I=1,IPARSIZ)

   10 FORMAT ('lda=',i6,/,'n=',i6,/,'blksz=',i6,/,'loclen=',i4,/,
     +       'basisdim=',i4,/,'nprocs=',i4,/,'procid=',i4,/,
     +       'precontype=',i4,/,'stoptype=',i4,/,'maxit=',i4,/,'itno=',
     +       i4,/,'status=',i4,/,'steperr=',i4)

      WRITE (6,FMT=20) (DPAR(I),I=1,DPARSIZ)

   20 FORMAT ('epsilon=',d20.12,/,'exitnorm=',d20.12,/,
     +       'eigenvalues region:',/,4(d20.12,1x))

      RETURN

      END
      SUBROUTINE PIMDQMR(X,B,WRK,IPAR,DPAR,MATVEC,TMATVEC,PRECONL,
     +                   PRECONR,PDSUM,PDNRM,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=2)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM
      EXTERNAL PDNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS,TMATVEC
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION ABSGAMMA2,ABSTAU02,DEN,EPSILON,EPSILON0,EXITNORM,
     +                 GAMMA,GAMMA0,KAPPA,KAPPA0,KSI,KSI0,LAMBDA,
     +                 LAMBDA0,MU,MU0,RHO,RHO0,RHSSTOP,TAU,TAU0,THETA
      INTEGER BASISDIM,BLKSZ,CNVRTX,IAP,IATWTILDE,ID,IP,IQ,IR,IS,ITNO,
     +        IVTILDE,IWTILDE,IXOLD,IZ,LDA,LOCLEN,MAXIT,N,NPROCS,
     +        PRECONTYPE,PROCID,STATUS,STEPERR,STOPTYPE
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DOTS(4)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT,DSETRHSSTOP
      EXTERNAL DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,DINIT,DSCAL,PIMDGETPAR,STOPCRIT
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS,SQRT
*     ..
      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 20

      END IF

*  Needs no conversion Y=Q2X for residual
      CNVRTX = 0

*  Set indices for mapping local vectors into wrk
      IR = 1
      IP = IR + LOCLEN
      IAP = IP + LOCLEN
      IQ = IAP + LOCLEN
      ID = IQ + LOCLEN
      IS = ID + LOCLEN
      IVTILDE = IS + LOCLEN
      IWTILDE = IVTILDE + LOCLEN
      IATWTILDE = IWTILDE + LOCLEN
      IXOLD = IATWTILDE + LOCLEN
      IZ = IXOLD + LOCLEN

*  Set RHS of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PDNRM)

*  1. lambda=1, kappa=-1, theta=-1
      LAMBDA = ONE
      KAPPA = -ONE
      THETA = -ONE

*  2. wtilde=vtilde=r=Q1(b-AQ2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL MATVEC(X,WRK(IAP),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IAP),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL MATVEC(X,WRK(IAP),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IAP),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL PRECONR(X,WRK(IAP),IPAR)
              CALL MATVEC(WRK(IAP),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL DCOPY(LOCLEN,B,1,WRK(IP),1)
              CALL PRECONR(X,WRK(IAP),IPAR)
              CALL MATVEC(WRK(IAP),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IP),1)
              CALL PRECONL(WRK(IP),WRK(IR),IPAR)
          END IF

      ELSE
*     r has been set to Q1b in the call to dsetrhsstop
          IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL MATVEC(X,WRK(IAP),IPAR)
              CALL PRECONL(WRK(IAP),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IAP),IPAR)
              CALL PRECONL(WRK(IAP),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)
          END IF

      END IF

      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IWTILDE),1)
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IVTILDE),1)

*  3. p=q=d=s=0
      CALL DINIT(LOCLEN,ZERO,WRK(IP),1)
      CALL DINIT(LOCLEN,ZERO,WRK(IQ),1)
      CALL DINIT(LOCLEN,ZERO,WRK(ID),1)
      CALL DINIT(LOCLEN,ZERO,WRK(IS),1)

*  4. gamma=||vtilde||_{2}, ksi=||wtilde||_{2},
*     rho=wtilde^{T}vtilde, epsilon=(A^{T}wtilde)^{T}vtilde, mu=0
      DOTS(1) = DDOT(LOCLEN,WRK(IVTILDE),1,WRK(IVTILDE),1)
      DOTS(2) = DDOT(LOCLEN,WRK(IWTILDE),1,WRK(IWTILDE),1)
      DOTS(3) = DDOT(LOCLEN,WRK(IWTILDE),1,WRK(IVTILDE),1)

*  Compute A^{T}wtilde
      IF (PRECONTYPE.EQ.0) THEN
          CALL TMATVEC(WRK(IWTILDE),WRK(IATWTILDE),IPAR)

      ELSE IF (PRECONTYPE.EQ.1) THEN
          CALL TMATVEC(WRK(IWTILDE),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IATWTILDE),IPAR)

      ELSE IF (PRECONTYPE.EQ.2) THEN
          CALL PRECONR(WRK(IWTILDE),WRK(IZ),IPAR)
          CALL TMATVEC(WRK(IZ),WRK(IATWTILDE),IPAR)

      ELSE IF (PRECONTYPE.EQ.3) THEN
          CALL PRECONR(WRK(IWTILDE),WRK(IATWTILDE),IPAR)
          CALL TMATVEC(WRK(IATWTILDE),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IATWTILDE),IPAR)
      END IF

      DOTS(4) = DDOT(LOCLEN,WRK(IATWTILDE),1,WRK(IVTILDE),1)

*  Accumulate simultaneously partial inner-products
      CALL PDSUM(4,DOTS,IPAR)

      GAMMA = SQRT(DOTS(1))
      KSI = SQRT(DOTS(2))
      RHO = DOTS(3)
      EPSILON = DOTS(4)
      MU = ZERO

*  5. tau=epsilon/rho
      IF (RHO.EQ.ZERO) THEN
          ITNO = 0
          STATUS = -3
          STEPERR = 5
          GO TO 20
      END IF

      TAU = EPSILON/RHO

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      DO 10 ITNO = 1,MAXIT

*  6. p=1/gamma*vtilde-mu*p
          IF (GAMMA.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 6
              GO TO 20

          END IF
          CALL DSCAL(LOCLEN,-MU,WRK(IP),1)
          CALL DAXPY(LOCLEN,ONE/GAMMA,WRK(IVTILDE),1,WRK(IP),1)

*  7. q=1/ksi*A^{T}wtilde-(gamma*mu)/ksi*q
          IF (KSI.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 7
              GO TO 20

          END IF
          CALL DSCAL(LOCLEN,- (GAMMA*MU)/KSI,WRK(IQ),1)
          CALL DAXPY(LOCLEN,ONE/KSI,WRK(IATWTILDE),1,WRK(IQ),1)

* 8. vtilde=Ap-tau/gamma*vtilde
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IP),WRK(IAP),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(WRK(IP),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IAP),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IP),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IAP),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IP),WRK(IAP),IPAR)
              CALL MATVEC(WRK(IAP),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IAP),IPAR)
          END IF

          CALL DSCAL(LOCLEN,-TAU/GAMMA,WRK(IVTILDE),1)
          CALL DAXPY(LOCLEN,ONE,WRK(IAP),1,WRK(IVTILDE),1)

*  9. wtilde=q-tau/ksi*wtilde
          CALL DSCAL(LOCLEN,-TAU/KSI,WRK(IWTILDE),1)
          CALL DAXPY(LOCLEN,ONE,WRK(IQ),1,WRK(IWTILDE),1)

* 10. check stopping criterion
          CALL DCOPY(LOCLEN,X,1,WRK(IXOLD),1)

          CALL STOPCRIT(B,WRK(IR),WRK(IZ),X,WRK(IXOLD),WRK(IATWTILDE),
     +                  RHSSTOP,CNVRTX,EXITNORM,STATUS,IPAR,MATVEC,
     +                  TMATVEC,PRECONR,PDSUM,PDNRM)

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IZ))

          IF (STATUS.EQ.-5) THEN
              STEPERR = 10
              GO TO 20
          ELSE IF (STATUS.EQ.0) THEN
              GO TO 20
          END IF

* 11. gamma=||vtilde||_{2}, ksi=||wtilde||_{2},
*     rho=wtilde^{T}vtilde, epsilon=(A^{T}wtilde)^{T}vtilde
          DOTS(1) = DDOT(LOCLEN,WRK(IVTILDE),1,WRK(IVTILDE),1)
          DOTS(2) = DDOT(LOCLEN,WRK(IWTILDE),1,WRK(IWTILDE),1)
          DOTS(3) = DDOT(LOCLEN,WRK(IWTILDE),1,WRK(IVTILDE),1)

*  Compute A^{T}wtilde
          IF (PRECONTYPE.EQ.0) THEN
              CALL TMATVEC(WRK(IWTILDE),WRK(IATWTILDE),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL TMATVEC(WRK(IWTILDE),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IATWTILDE),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IWTILDE),WRK(IZ),IPAR)
              CALL TMATVEC(WRK(IZ),WRK(IATWTILDE),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IWTILDE),WRK(IATWTILDE),IPAR)
              CALL TMATVEC(WRK(IATWTILDE),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IATWTILDE),IPAR)
          END IF

          DOTS(4) = DDOT(LOCLEN,WRK(IATWTILDE),1,WRK(IVTILDE),1)

*  Accumulate simultaneously partial inner-products
          CALL PDSUM(4,DOTS,IPAR)

          GAMMA0 = GAMMA
          GAMMA = SQRT(DOTS(1))
          KSI0 = KSI
          KSI = SQRT(DOTS(2))
          RHO0 = RHO
          RHO = DOTS(3)
          EPSILON0 = EPSILON
          EPSILON = DOTS(4)

* 12. mu=(gamma0*ksi0*rho)/(gamma*tau*rho0)
          DEN = GAMMA*TAU*RHO0
          IF (DEN.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 12
              GO TO 20

          END IF
          MU0 = MU
          MU = (GAMMA0*KSI0*RHO)/DEN

* 13. tau=epsilon/rho-gamma*mu
          IF (RHO.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 13
              GO TO 20

          END IF
          TAU0 = TAU
          TAU = EPSILON/RHO - GAMMA*MU

* 14. theta=(|tau0|^2*(1-lambda))/(lambda*|tau|^2+|gamma|^2)
          ABSTAU02 = ABS(TAU0)**2
          ABSGAMMA2 = ABS(GAMMA)**2
          DEN = LAMBDA*ABSTAU02 + ABSGAMMA2
          IF (DEN.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 14
              GO TO 20

          END IF
          THETA = (ABSTAU02* (ONE-LAMBDA))/DEN

* 15. kappa=(-gamma0*CONJG(tau0)*kappa0)/(gamma0*|tau|^2+|gamma|^2)
          KAPPA0 = KAPPA
          KAPPA = - (GAMMA0*TAU0*KAPPA0)/DEN

* 16. lambda=(lambda0*|tau0|^2)/(gamma0*|tau|^2+|gamma|^2)
          LAMBDA0 = LAMBDA
          LAMBDA = LAMBDA0*ABSTAU02/DEN

* 17. d=theta*d+kappa*p
          CALL DSCAL(LOCLEN,THETA,WRK(ID),1)
          CALL DAXPY(LOCLEN,KAPPA,WRK(IP),1,WRK(ID),1)

* 18. s=theta*s+kappa*A*p
          CALL DSCAL(LOCLEN,THETA,WRK(IS),1)
          CALL DAXPY(LOCLEN,KAPPA,WRK(IAP),1,WRK(IS),1)

* 19. x=x+d
          CALL DAXPY(LOCLEN,ONE,WRK(ID),1,X,1)

* 20. r=r-s
          CALL DAXPY(LOCLEN,-ONE,WRK(IS),1,WRK(IR),1)

   10 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

   20 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMDRBICGSTAB(X,B,WRK,IPAR,DPAR,MATVEC,PRECONL,PRECONR,
     +                         PDSUM,PDNRM2,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=2)
      INTEGER IBDIM
      PARAMETER (IBDIM=10)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM2
      EXTERNAL PDNRM2
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION ALPHA,BETA,EPSILON,EXITNORM,KSI,OMEGA,RHO0,RHO1,
     +                 RHSSTOP,S
      INTEGER BASISDIM,BLKSZ,CNVRTX,I,I0,I1,I2,I3,I4,IR,IRTILDE,ITNO,IU,
     +        IW,IXOLD,IZ,J,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,PROCID,
     +        STATUS,STEPERR,STOPTYPE
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DOTS(IBDIM),GAMMA(IBDIM),GAMMA1(IBDIM),
     +                 GAMMA2(IBDIM),SIGMA(IBDIM),TAU(IBDIM,IBDIM)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT,DSETRHSSTOP
      EXTERNAL DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,DINIT,PIMDGETPAR,STOPCRIT
*     ..
      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 9999

      END IF

*  Does not need conversion Y=Q2X for residual
      CNVRTX = 0

*  Set indices for mapping local vectors into wrk
      IRTILDE = 1
      IW = IRTILDE + LOCLEN
      IZ = IW + LOCLEN
      IXOLD = IZ + LOCLEN
      IR = IXOLD + LOCLEN
      IU = IR + (BASISDIM+1)*LOCLEN

*  Set rhs of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PDNRM2)

*  1. r=Q1(b-AQ2x)
      IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
          CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
          CALL MATVEC(X,WRK(IW),IPAR)
          CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IR),1)

      ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
          CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
          CALL MATVEC(X,WRK(IW),IPAR)
          CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
          CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

      ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
          CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
          CALL PRECONR(X,WRK(IW),IPAR)
          CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
          CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

      ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
          CALL DCOPY(LOCLEN,B,1,WRK(IW),1)
          CALL PRECONR(X,WRK(IR),IPAR)
          CALL MATVEC(WRK(IR),WRK(IZ),IPAR)
          CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IW),1)
          CALL PRECONL(WRK(IW),WRK(IR),IPAR)
      END IF

*  2. rtilde=r
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IRTILDE),1)

*  3. u0=0
      CALL DINIT(LOCLEN,ZERO,WRK(IU),1)

*  4. rho0=1, alpha=0, omega=1
      RHO0 = ONE
      ALPHA = ZERO
      OMEGA = ONE

*  Loop
      STATUS = 0
      STEPERR = -1
      EXITNORM = -ONE
      DO 120 ITNO = 1,MAXIT

*  5. rho0=-omega*rho0
          RHO0 = -OMEGA*RHO0

*  BiCG loop
          I1 = 0
          I2 = LOCLEN
          DO 30 J = 0,BASISDIM - 1

*  6. rho1=r(j)^{T}rtilde
              DOTS(1) = DDOT(LOCLEN,WRK(IR+I1),1,WRK(IRTILDE),1)
              CALL PDSUM(1,DOTS,IPAR)
              RHO1 = DOTS(1)

*  7. beta=alpha*rho1/rho0
              IF (RHO0.EQ.ZERO) THEN
                  STATUS = -3
                  STEPERR = 7
                  GO TO 9999

              END IF

              BETA = ALPHA*RHO1/RHO0

*  8. rho0=rho1
              RHO0 = RHO1

*  9. u(i)=r(i)-beta*u(i), i=0:j
              I3 = 0
              DO 10 I = 0,J
                  CALL DCOPY(LOCLEN,WRK(IU+I3),1,WRK(IZ),1)
                  CALL DCOPY(LOCLEN,WRK(IR+I3),1,WRK(IU+I3),1)
                  CALL DAXPY(LOCLEN,-BETA,WRK(IZ),1,WRK(IU+I3),1)
                  I3 = I3 + LOCLEN
   10         CONTINUE

* 10. u(j+1)=Q1AQ2u(j)
              IF (PRECONTYPE.EQ.0) THEN
                  CALL MATVEC(WRK(IU+I1),WRK(IU+I2),IPAR)

              ELSE IF (PRECONTYPE.EQ.1) THEN
                  CALL MATVEC(WRK(IU+I1),WRK(IW),IPAR)
                  CALL PRECONL(WRK(IW),WRK(IU+I2),IPAR)

              ELSE IF (PRECONTYPE.EQ.2) THEN
                  CALL PRECONR(WRK(IU+I1),WRK(IW),IPAR)
                  CALL MATVEC(WRK(IW),WRK(IU+I2),IPAR)

              ELSE IF (PRECONTYPE.EQ.3) THEN
                  CALL PRECONR(WRK(IU+I1),WRK(IZ),IPAR)
                  CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
                  CALL PRECONL(WRK(IW),WRK(IU+I2),IPAR)
              END IF

* 11. ksi=u(j+1)^{T}rtilde
              DOTS(1) = DDOT(LOCLEN,WRK(IU+I2),1,WRK(IRTILDE),1)
              CALL PDSUM(1,DOTS,IPAR)
              KSI = DOTS(1)

* 12. alpha=rho0/ksi
              IF (KSI.EQ.ZERO) THEN
                  STATUS = -3
                  STEPERR = 12
                  GO TO 9999

              END IF

              ALPHA = RHO0/KSI

* 13. r(i)=r(i)-alpha*u(i+1), i=0:j
              I3 = 0
              I4 = LOCLEN
              DO 20 I = 0,J
                  CALL DAXPY(LOCLEN,-ALPHA,WRK(IU+I4),1,WRK(IR+I3),1)
                  I3 = I3 + LOCLEN
                  I4 = I4 + LOCLEN
   20         CONTINUE

* 14. r(j+1)=Q1AQ2r(j)
              IF (PRECONTYPE.EQ.0) THEN
                  CALL MATVEC(WRK(IR+I1),WRK(IR+I2),IPAR)

              ELSE IF (PRECONTYPE.EQ.1) THEN
                  CALL MATVEC(WRK(IR+I1),WRK(IW),IPAR)
                  CALL PRECONL(WRK(IW),WRK(IR+I2),IPAR)

              ELSE IF (PRECONTYPE.EQ.2) THEN
                  CALL PRECONR(WRK(IR+I1),WRK(IW),IPAR)
                  CALL MATVEC(WRK(IW),WRK(IR+I2),IPAR)

              ELSE IF (PRECONTYPE.EQ.3) THEN
                  CALL PRECONR(WRK(IR+I1),WRK(IZ),IPAR)
                  CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
                  CALL PRECONL(WRK(IW),WRK(IR+I2),IPAR)
              END IF

* 15. x0=x0+alpha*u0
              CALL DCOPY(LOCLEN,X,1,WRK(IXOLD),1)
              CALL DAXPY(LOCLEN,ALPHA,WRK(IU),1,X,1)
              I1 = I1 + LOCLEN
              I2 = I2 + LOCLEN
   30     CONTINUE

* 16. check stopping criterion
          CALL STOPCRIT(B,WRK(IR),WRK(IZ),X,WRK(IXOLD),WRK(IW),RHSSTOP,
     +                  CNVRTX,EXITNORM,STATUS,IPAR,MATVEC,MATVEC,
     +                  PRECONR,PDSUM,PDNRM2)

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IZ))

          IF (STATUS.EQ.-5) THEN
              STEPERR = 16
              GO TO 9999
          ELSE IF (STATUS.EQ.0) THEN
              GO TO 9999
          END IF
*  MR loop

* 17. sigma(1)=r(1)^{T}r(1), gamma'(1)=r(0)^{T}r(1)/sigma(1)
          DOTS(1) = DDOT(LOCLEN,WRK(IR+LOCLEN),1,WRK(IR+LOCLEN),1)
          DOTS(2) = DDOT(LOCLEN,WRK(IR),1,WRK(IR+LOCLEN),1)
          CALL PDSUM(2,DOTS,IPAR)
          SIGMA(1) = DOTS(1)

          IF (SIGMA(1).EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 17
              GO TO 9999

          END IF

          GAMMA1(1) = DOTS(2)/SIGMA(1)

          I0 = LOCLEN + LOCLEN
          DO 60 J = 2,BASISDIM

* 18. tau(i,j)=r(j)^{T}r(i)/sigma(i), r(j)=r(j)-tau(i,j)r(i)
              I1 = LOCLEN
              DO 40 I = 1,J - 1
                  DOTS(I) = DDOT(LOCLEN,WRK(IR+I0),1,WRK(IR+I1),1)
                  I1 = I1 + LOCLEN
   40         CONTINUE
              CALL PDSUM(J-1,DOTS,IPAR)
              I1 = LOCLEN
              DO 50 I = 1,J - 1
                  TAU(I,J) = DOTS(I)/SIGMA(I)
                  CALL DAXPY(LOCLEN,-TAU(I,J),WRK(IR+I1),1,WRK(IR+I0),1)
   50         CONTINUE

* 19. sigma(j)=r(j)^{T}r(j), gamma'(j)=r(0)^{T}r(j)/sigma(j)
              DOTS(1) = DDOT(LOCLEN,WRK(IR+I0),1,WRK(IR+I0),1)
              DOTS(2) = DDOT(LOCLEN,WRK(IR),1,WRK(IR+I0),1)
              CALL PDSUM(2,DOTS,IPAR)
              SIGMA(J) = DOTS(1)

              IF (SIGMA(J).EQ.ZERO) THEN
                  STATUS = -3
                  STEPERR = 19
                  GO TO 9999

              END IF

              GAMMA1(J) = DOTS(2)/SIGMA(J)
              I0 = I0 + LOCLEN
   60     CONTINUE

* 20. gamma_{l}=omega=gamma'_{l}
*     gamma_{j}=gamma'_{j}-\sum_{i=j+1}^{l}{tau_{j,i}gamma_{i}}
          GAMMA(BASISDIM) = GAMMA1(BASISDIM)
          OMEGA = GAMMA(BASISDIM)
          DO 80 J = BASISDIM - 1,1,-1
              S = ZERO
              DO 70 I = J + 1,BASISDIM
                  S = S + TAU(J,I)*GAMMA(I)
   70         CONTINUE
              GAMMA(J) = GAMMA1(J) - S
   80     CONTINUE

* 21. gamma''=gamma_{j+1}+\sum_{i=j+1}^{l-1}{tau_{j,i}gamma_{i+1}}
          DO 100 J = 1,BASISDIM - 1
              S = ZERO
              DO 90 I = J + 1,BASISDIM - 1
                  S = S + TAU(J,I)*GAMMA(I+1)
   90         CONTINUE
              GAMMA2(J) = GAMMA(J+1) + S
  100     CONTINUE

*  Update

* 22. x(0)=x(0)+gamma(1)r(0)
          CALL DAXPY(LOCLEN,GAMMA(1),WRK(IR),1,X,1)

* 23. r(0)=r(0)-gamma'(l)r(l)
          CALL DAXPY(LOCLEN,-GAMMA1(BASISDIM),WRK(IR+BASISDIM*LOCLEN),1,
     +               WRK(IR),1)

* 24. u(0)=u(0)-gamma(l)u(l)
          CALL DAXPY(LOCLEN,-GAMMA(BASISDIM),WRK(IU+BASISDIM*LOCLEN),1,
     +               WRK(IU),1)

          I0 = LOCLEN
          DO 110 J = 1,BASISDIM - 1

* 25. u(0)=u(0)-gamma(j)u(j), j=1:l-1
              CALL DAXPY(LOCLEN,-GAMMA(J),WRK(IU+I0),1,WRK(IU),1)

* 26. x(0)=x(0)+gamma''(j)r(j), j=1:l-1
              CALL DAXPY(LOCLEN,GAMMA2(J),WRK(IR+I0),1,X,1)

* 27. r(0)=r(0)-gamma'(j)r(j), j=1:l-1
              CALL DAXPY(LOCLEN,-GAMMA1(J),WRK(IR+I0),1,WRK(IR),1)
              I0 = I0 + LOCLEN
  110     CONTINUE

  120 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMDRGCR(X,B,WRK,IPAR,DPAR,MATVEC,PRECONL,PRECONR,
     +                    PDSUM,PDNRM,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=2)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM
      EXTERNAL PDNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION ALPHA,BETA,EPSILON,EXITNORM,RHSSTOP,XI
      INTEGER BASISDIM,BLKSZ,CNVRTX,I,IDOTS,IP,IQ,IR,ITNO,IW,IXOLD,IZ,
     +        IZETA,J,J0,K,K1,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,
     +        PROCID,STATUS,STEPERR,STOPTYPE
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DOTS(2)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT,DSETRHSSTOP
      EXTERNAL DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,DINIT,PIMDGETPAR,STOPCRIT
*     ..
      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 9999

      END IF

*  Does not need conversion Y=Q2X for residual
      CNVRTX = 0

*  Set indices for mapping local vectors into wrk
      IR = 1
      IP = IR + LOCLEN
      IW = IP + BASISDIM*LOCLEN
      IZETA = IW + BASISDIM*LOCLEN
      IZ = IZETA + BASISDIM
      IQ = IZ + LOCLEN
      IXOLD = IQ + LOCLEN
      IDOTS = IXOLD + LOCLEN

*  Set rhs of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PDNRM)

*  1. r=Q1(b-AQ2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL DCOPY(LOCLEN,B,1,WRK(IP),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IP),1)
              CALL PRECONL(WRK(IP),WRK(IR),IPAR)
          END IF

      ELSE
*     r has been set to Qb in the call to dsetrhsstop
          IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)
          END IF

      END IF

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      DO 20 ITNO = 1,MAXIT

*  2. p(1)=r
          CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP),1)

          K = 0
          DO 40 J = 1,BASISDIM
              J0 = J - 1

*  3. w(j)=Q1AQ2p(j)
              IF (PRECONTYPE.EQ.0) THEN
                  CALL MATVEC(WRK(IP+K),WRK(IW+K),IPAR)

              ELSE IF (PRECONTYPE.EQ.1) THEN
                  CALL MATVEC(WRK(IP+K),WRK(IZ),IPAR)
                  CALL PRECONL(WRK(IZ),WRK(IW+K),IPAR)

              ELSE IF (PRECONTYPE.EQ.2) THEN
                  CALL PRECONR(WRK(IP+K),WRK(IZ),IPAR)
                  CALL MATVEC(WRK(IZ),WRK(IW+K),IPAR)

              ELSE IF (PRECONTYPE.EQ.3) THEN
                  CALL PRECONR(WRK(IP+K),WRK(IW+K),IPAR)
                  CALL MATVEC(WRK(IW+K),WRK(IZ),IPAR)
                  CALL PRECONL(WRK(IZ),WRK(IW+K),IPAR)
              END IF

*  4. zeta(j)=dot(w(j),w(j))
              DOTS(1) = DDOT(LOCLEN,WRK(IW+K),1,WRK(IW+K),1)

*  5. xi=dot(r,w(j))
              DOTS(2) = DDOT(LOCLEN,WRK(IR),1,WRK(IW+K),1)

*  Accumulate simultaneously partial values
              CALL PDSUM(2,DOTS,IPAR)
              WRK(IZETA+J0) = DOTS(1)
              XI = DOTS(2)

*  5. alpha=xi/zeta(j)
              IF (WRK(IZETA+J0).EQ.ZERO) THEN
                  STATUS = -3
                  STEPERR = 5
                  GO TO 9999

              END IF

              ALPHA = XI/WRK(IZETA+J0)

*  6. x=x+alpha*p(j)
              CALL DCOPY(LOCLEN,X,1,WRK(IXOLD),1)
              CALL DAXPY(LOCLEN,ALPHA,WRK(IP+K),1,X,1)

*  7. r=r-alpha*w(j)
              CALL DAXPY(LOCLEN,-ALPHA,WRK(IW+K),1,WRK(IR),1)

*  8. check stopping criterion
              CALL STOPCRIT(B,WRK(IR),WRK(IZ),X,WRK(IXOLD),WRK(IQ),
     +                      RHSSTOP,CNVRTX,EXITNORM,STATUS,IPAR,MATVEC,
     +                      MATVEC,PRECONR,PDSUM,PDNRM)

*  Call monitoring routine
              CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IZ))

              IF (STATUS.EQ.-5) THEN
                  STEPERR = 8
                  GO TO 9999
              ELSE IF (STATUS.EQ.0) THEN
                  GO TO 9999
              END IF
*  9. q=Q1AQ2r
              IF (PRECONTYPE.EQ.0) THEN
                  CALL MATVEC(WRK(IR),WRK(IQ),IPAR)

              ELSE IF (PRECONTYPE.EQ.1) THEN
                  CALL MATVEC(WRK(IR),WRK(IZ),IPAR)
                  CALL PRECONL(WRK(IZ),WRK(IQ),IPAR)

              ELSE IF (PRECONTYPE.EQ.2) THEN
                  CALL PRECONR(WRK(IR),WRK(IZ),IPAR)
                  CALL MATVEC(WRK(IZ),WRK(IQ),IPAR)

              ELSE IF (PRECONTYPE.EQ.3) THEN
                  CALL PRECONR(WRK(IR),WRK(IQ),IPAR)
                  CALL MATVEC(WRK(IQ),WRK(IZ),IPAR)
                  CALL PRECONL(WRK(IZ),WRK(IQ),IPAR)
              END IF

* 10. p(j+1)=r-sum_{i=1}^{j}{dot(q,w(i))/zeta(i)*p(i)}
              IF (J.LT.BASISDIM) THEN

*  Compute partial inner-products
                  K1 = 0
                  DO 50 I = 0,J - 1
                      WRK(IDOTS+I) = DDOT(LOCLEN,WRK(IQ),1,WRK(IW+K1),1)
                      K1 = K1 + LOCLEN
   50             CONTINUE

*  Accumulate simultaneously partial values
                  CALL PDSUM(J,WRK(IDOTS),IPAR)

*  Compute summation
                  CALL DINIT(LOCLEN,ZERO,WRK(IZ),1)
                  K1 = 0
                  DO 60 I = 0,J - 1
                      BETA = WRK(IDOTS+I)/WRK(IZETA+I)
                      CALL DAXPY(LOCLEN,BETA,WRK(IP+K1),1,WRK(IZ),1)
                      K1 = K1 + LOCLEN
   60             CONTINUE

*  Compute p(j+1)
                  K = K + LOCLEN
                  CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IP+K),1)
                  CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IP+K),1)
              END IF

   40     CONTINUE

   20 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMDRGMRES(X,B,WRK,IPAR,DPAR,MATVEC,PRECONL,PRECONR,
     +                      PDSUM,PDNRM2,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=2)
      INTEGER IBDIM
      PARAMETER (IBDIM=50)
      INTEGER LDR
      PARAMETER (LDR=IBDIM+1)
      INTEGER LDG
      PARAMETER (LDG=IBDIM+2)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM2
      EXTERNAL PDNRM2
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION BETA,EPSILON,ETA,EXITNORM,KSI,RHSSTOP,TAU1,TAU2
      INTEGER BASISDIM,BLKSZ,I,IRES,ITNO,IV,IW,IZ,J,K0,K1,LDA,LOCLEN,
     +        MAXIT,N,NPROCS,PRECONTYPE,PROCID,STATUS,STEPERR,STOPTYPE
      LOGICAL ENDED
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION G(LDG),R(LDR,LDR),RHO(IBDIM)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT,DSETRHSSTOP
      EXTERNAL DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,DECODE,DINIT,DSCAL,DTRSV,ENCODE,GIVENS,
     +         PIMDGETPAR
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS
*     ..
      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 9999

      END IF

*  Set indices for mapping local vectors into wrk
      IRES = 1
      IZ = IRES + LOCLEN
      IW = IZ + LOCLEN
      IV = IW + LOCLEN

*  Set rhs of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IRES),EPSILON,IPAR,PRECONL,PDNRM2)

*  1. r=Q1(b-AQ2x)
      IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
          CALL DCOPY(LOCLEN,B,1,WRK(IRES),1)
          CALL MATVEC(X,WRK(IW),IPAR)
          CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IRES),1)

      ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
          CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
          CALL MATVEC(X,WRK(IW),IPAR)
          CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
          CALL PRECONL(WRK(IZ),WRK(IRES),IPAR)

      ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
          CALL DCOPY(LOCLEN,B,1,WRK(IRES),1)
          CALL PRECONR(X,WRK(IW),IPAR)
          CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
          CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IRES),1)

      ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
          CALL DCOPY(LOCLEN,B,1,WRK(IW),1)
          CALL PRECONR(X,WRK(IRES),IPAR)
          CALL MATVEC(WRK(IRES),WRK(IZ),IPAR)
          CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IW),1)
          CALL PRECONL(WRK(IW),WRK(IRES),IPAR)
      END IF

*  2. beta=||r||_2
      BETA = PDNRM2(LOCLEN,WRK(IRES),IPAR)

*  Loop
      STATUS = 0
      STEPERR = -1
      EXITNORM = -ONE
      ENDED = .FALSE.
      DO 20 ITNO = 1,MAXIT

*  3. g=(beta,beta,...)
          G(1) = BETA
          G(2) = BETA

*  4. V(1)=r/beta
          IF (BETA.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 4
              GO TO 9999

          END IF

          CALL DCOPY(LOCLEN,WRK(IRES),1,WRK(IV),1)
          CALL DSCAL(LOCLEN,ONE/BETA,WRK(IV),1)

          K0 = 0
          DO 40 J = 1,BASISDIM

*     z=Q1AQ2V(j)
              IF (PRECONTYPE.EQ.0) THEN
                  CALL MATVEC(WRK(IV+K0),WRK(IZ),IPAR)

              ELSE IF (PRECONTYPE.EQ.1) THEN
                  CALL MATVEC(WRK(IV+K0),WRK(IW),IPAR)
                  CALL PRECONL(WRK(IW),WRK(IZ),IPAR)

              ELSE IF (PRECONTYPE.EQ.2) THEN
                  CALL PRECONR(WRK(IV+K0),WRK(IW),IPAR)
                  CALL MATVEC(WRK(IW),WRK(IZ),IPAR)

              ELSE IF (PRECONTYPE.EQ.3) THEN
                  CALL PRECONR(WRK(IV+K0),WRK(IZ),IPAR)
                  CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
                  CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              END IF


*  5. R(i,j)=dot(V(i),Q1AQ2V(j))
              K1 = 0
              DO 50 I = 1,J
                  R(I,J) = DDOT(LOCLEN,WRK(IV+K1),1,WRK(IZ),1)
                  K1 = K1 + LOCLEN
   50         CONTINUE
              CALL PDSUM(J,R(1,J),IPAR)

*  6. Vhat(j)=Q1AQ2V(j)-sum_{i=1}^{j}{R(i,j)V(i)}
              K1 = 0
              CALL DINIT(LOCLEN,ZERO,WRK(IW),1)
              DO 60 I = 1,J
                  CALL DAXPY(LOCLEN,R(I,J),WRK(IV+K1),1,WRK(IW),1)
                  K1 = K1 + LOCLEN
   60         CONTINUE
              CALL DSCAL(LOCLEN,-ONE,WRK(IW),1)
              CALL DAXPY(LOCLEN,ONE,WRK(IZ),1,WRK(IW),1)

*  From this point, w holds the (j+1)-st column of vhat

*  7. R(j+1,j)=||Vhat(j)||_2
              R(J+1,J) = PDNRM2(LOCLEN,WRK(IW),IPAR)

*  8. V(j+1)=Vhat(j)/R(j+1,j)
              IF (R(J+1,J).EQ.ZERO) THEN
                  STATUS = -2
                  STEPERR = 8
                  GO TO 9999

              END IF

              K0 = K0 + LOCLEN
              CALL DSCAL(LOCLEN,ONE/R(J+1,J),WRK(IW),1)
              CALL DCOPY(LOCLEN,WRK(IW),1,WRK(IV+K0),1)

*  9. Apply previous Givens' rotations to column j of R
              DO 70 I = 1,J - 1
                  CALL DECODE(RHO(I),KSI,ETA)
                  TAU1 = R(I,J)
                  TAU2 = R(I+1,J)
                  R(I,J) = KSI*TAU1 - ETA*TAU2
                  R(I+1,J) = ETA*TAU1 + KSI*TAU2
   70         CONTINUE

* 10. Compute Givens' rotation to zero element R(j+1,j)
              CALL GIVENS(R(J,J),R(J+1,J),KSI,ETA)
              TAU1 = R(J,J)
              TAU2 = R(J+1,J)
              R(J,J) = KSI*TAU1 - ETA*TAU2
              R(J+1,J) = ETA*TAU1 + KSI*TAU2
              CALL ENCODE(RHO(J),KSI,ETA)

*  11. Update g
              G(J) = G(J)*KSI
              G(J+1) = G(J+1)*ETA
              G(J+2) = G(J+1)

*  12. If |g(j+1)|<rhsstop stop
              EXITNORM = ABS(G(J+1))
              IF (EXITNORM.LT.RHSSTOP) THEN
                  BASISDIM = J
                  ENDED = .TRUE.
                  GO TO 80

              END IF

   40     CONTINUE
   80     CONTINUE

*  13. Solve Ry=g (solution to least-squares problem)
          CALL DTRSV('U','N','N',BASISDIM,R,LDR,G,1)

*  14. x=x+Vy (Form approximate solution after a c-cycle)
          K1 = 0
          DO 100 I = 1,BASISDIM
              CALL DAXPY(LOCLEN,G(I),WRK(IV+K1),1,X,1)
              K1 = K1 + LOCLEN
  100     CONTINUE

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IRES),WRK(IRES))

          IF (ENDED) GO TO 9999

*  15. r=Q1(b-AQ2x)
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL DCOPY(LOCLEN,B,1,WRK(IRES),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IRES),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IRES),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
              CALL DCOPY(LOCLEN,B,1,WRK(IRES),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IRES),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL DCOPY(LOCLEN,B,1,WRK(IW),1)
              CALL PRECONR(X,WRK(IRES),IPAR)
              CALL MATVEC(WRK(IRES),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IW),1)
              CALL PRECONL(WRK(IW),WRK(IRES),IPAR)
          END IF

*  16. beta=||r||_2
          BETA = PDNRM2(LOCLEN,WRK(IRES),IPAR)

   20 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMDRGMRESEV(X,B,WRK,IPAR,DPAR,MATVEC,PRECONL,PRECONR,
     +                        PDSUM,PDNRM2,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=6)
      INTEGER IBDIM
      PARAMETER (IBDIM=50)
      INTEGER LDR
      PARAMETER (LDR=IBDIM+1)
      INTEGER LDG
      PARAMETER (LDG=IBDIM+2)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM2
      EXTERNAL PDNRM2
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION BETA,EPSILON,ETA,EXITNORM,KSI,RHSSTOP,TAU1,TAU2
      INTEGER BASISDIM,BLKSZ,I,INFO,IRES,ITNO,IV,IW,IZ,J,K0,K1,LDA,
     +        LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,PROCID,STATUS,STEPERR,
     +        STOPTYPE
      LOGICAL ENDED
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DUMMY(1,1),EIGHI(IBDIM),EIGHR(IBDIM),G(LDG),
     +                 H(LDR,IBDIM),R(LDR,IBDIM),RHO(IBDIM),WRKEV(IBDIM)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT,DSETRHSSTOP
      EXTERNAL DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,DECODE,DHSEQR,DINIT,DSCAL,DTRSV,ENCODE,
     +         GIVENS,PIMDGETPAR
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS,MAX,MIN
*     ..
      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 9999

      END IF

*  Set indices for mapping local vectors into wrk
      IRES = 1
      IZ = IRES + LOCLEN
      IW = IZ + LOCLEN
      IV = IW + LOCLEN

*  Set rhs of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IRES),EPSILON,IPAR,PRECONL,PDNRM2)

*  1. r=Q1(b-AQ2x)
      IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
          CALL DCOPY(LOCLEN,B,1,WRK(IRES),1)
          CALL MATVEC(X,WRK(IW),IPAR)
          CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IRES),1)

      ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
          CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
          CALL MATVEC(X,WRK(IW),IPAR)
          CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
          CALL PRECONL(WRK(IZ),WRK(IRES),IPAR)

      ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
          CALL DCOPY(LOCLEN,B,1,WRK(IRES),1)
          CALL PRECONR(X,WRK(IW),IPAR)
          CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
          CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IRES),1)

      ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
          CALL DCOPY(LOCLEN,B,1,WRK(IW),1)
          CALL PRECONR(X,WRK(IRES),IPAR)
          CALL MATVEC(WRK(IRES),WRK(IZ),IPAR)
          CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IW),1)
          CALL PRECONL(WRK(IW),WRK(IRES),IPAR)
      END IF

*  2. beta=||r||_2
      BETA = PDNRM2(LOCLEN,WRK(IRES),IPAR)

*  Loop
      STATUS = 0
      STEPERR = -1
      EXITNORM = -ONE
      ENDED = .FALSE.
      DO 20 ITNO = 1,MAXIT

*  3. g=(beta,beta,...)
          G(1) = BETA
          G(2) = BETA

*  4. V(1)=r/beta
          IF (BETA.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 4
              GO TO 9999

          END IF

          CALL DCOPY(LOCLEN,WRK(IRES),1,WRK(IV),1)
          CALL DSCAL(LOCLEN,ONE/BETA,WRK(IV),1)

          K0 = 0
          DO 40 J = 1,BASISDIM

*     z=Q1AQ2V(j)
              IF (PRECONTYPE.EQ.0) THEN
                  CALL MATVEC(WRK(IV+K0),WRK(IZ),IPAR)

              ELSE IF (PRECONTYPE.EQ.1) THEN
                  CALL MATVEC(WRK(IV+K0),WRK(IW),IPAR)
                  CALL PRECONL(WRK(IW),WRK(IZ),IPAR)

              ELSE IF (PRECONTYPE.EQ.2) THEN
                  CALL PRECONR(WRK(IV+K0),WRK(IW),IPAR)
                  CALL MATVEC(WRK(IW),WRK(IZ),IPAR)

              ELSE IF (PRECONTYPE.EQ.3) THEN
                  CALL PRECONR(WRK(IV+K0),WRK(IZ),IPAR)
                  CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
                  CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              END IF


*  5. R(i,j)=dot(V(i),Q1AQ2V(j))
              K1 = 0
              DO 50 I = 1,J
                  H(I,J) = DDOT(LOCLEN,WRK(IV+K1),1,WRK(IZ),1)
                  R(I,J) = H(I,J)
                  K1 = K1 + LOCLEN
   50         CONTINUE
              CALL PDSUM(J,R(1,J),IPAR)

*  6. Vhat(j)=Q1AQ2V(j)-sum_{i=1}^{j}{R(i,j)V(i)}
              K1 = 0
              CALL DINIT(LOCLEN,ZERO,WRK(IW),1)
              DO 60 I = 1,J
                  CALL DAXPY(LOCLEN,R(I,J),WRK(IV+K1),1,WRK(IW),1)
                  K1 = K1 + LOCLEN
   60         CONTINUE
              CALL DSCAL(LOCLEN,-ONE,WRK(IW),1)
              CALL DAXPY(LOCLEN,ONE,WRK(IZ),1,WRK(IW),1)

*  From this point, w holds the (j+1)-st column of vhat

*  7. R(j+1,j)=||Vhat(j)||_2
              H(J+1,J) = PDNRM2(LOCLEN,WRK(IW),IPAR)
              R(J+1,J) = H(J+1,J)

*  8. V(j+1)=Vhat(j)/R(j+1,j)
              IF (R(J+1,J).EQ.ZERO) THEN
                  STATUS = -2
                  STEPERR = 8
                  GO TO 9999

              END IF

              K0 = K0 + LOCLEN
              CALL DSCAL(LOCLEN,ONE/R(J+1,J),WRK(IW),1)
              CALL DCOPY(LOCLEN,WRK(IW),1,WRK(IV+K0),1)

*  9. Apply previous Givens' rotations to column j of R
              DO 70 I = 1,J - 1
                  CALL DECODE(RHO(I),KSI,ETA)
                  TAU1 = R(I,J)
                  TAU2 = R(I+1,J)
                  R(I,J) = KSI*TAU1 - ETA*TAU2
                  R(I+1,J) = ETA*TAU1 + KSI*TAU2
   70         CONTINUE

* 10. Compute Givens' rotation to zero element R(j+1,j)
              CALL GIVENS(R(J,J),R(J+1,J),KSI,ETA)
              TAU1 = R(J,J)
              TAU2 = R(J+1,J)
              R(J,J) = KSI*TAU1 - ETA*TAU2
              R(J+1,J) = ETA*TAU1 + KSI*TAU2
              CALL ENCODE(RHO(J),KSI,ETA)

*  11. Update g
              G(J) = G(J)*KSI
              G(J+1) = G(J+1)*ETA
              G(J+2) = G(J+1)

*  12. If |g(j+1)|<rhsstop stop
              EXITNORM = ABS(G(J+1))
              IF (EXITNORM.LT.RHSSTOP) THEN
                  BASISDIM = J
                  ENDED = .TRUE.
                  GO TO 80

              END IF

   40     CONTINUE
   80     CONTINUE

*  13. Solve Ry=g (solution to least-squares problem)
          CALL DTRSV('U','N','N',BASISDIM,R,LDR,G,1)

*  14. x=x+Vy (Form approximate solution after a c-cycle)
          K1 = 0
          DO 100 I = 1,BASISDIM
              CALL DAXPY(LOCLEN,G(I),WRK(IV+K1),1,X,1)
              K1 = K1 + LOCLEN
  100     CONTINUE

*  15. compute eigenvalues of H, obtain bounding region
          CALL DHSEQR('E','N',BASISDIM,1,BASISDIM,H,LDR,EIGHR,EIGHI,
     +                DUMMY,1,WRKEV,IBDIM,INFO)
          DPAR(3) = EIGHR(1)
          DPAR(4) = EIGHR(1)
          DPAR(5) = EIGHI(1)
          DPAR(6) = EIGHI(1)
          DO 110 I = 2,BASISDIM
              DPAR(3) = MIN(DPAR(3),EIGHR(I))
              DPAR(4) = MAX(DPAR(4),EIGHR(I))
              DPAR(5) = MIN(DPAR(5),EIGHI(I))
              DPAR(6) = MAX(DPAR(6),EIGHI(I))
  110     CONTINUE

*  Call monitoring routine
          CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IRES),WRK(IRES))

          IF (ENDED) GO TO 9999

*  16. r=Q1(b-AQ2x)
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL DCOPY(LOCLEN,B,1,WRK(IRES),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IRES),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IRES),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
              CALL DCOPY(LOCLEN,B,1,WRK(IRES),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IRES),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL DCOPY(LOCLEN,B,1,WRK(IW),1)
              CALL PRECONR(X,WRK(IRES),IPAR)
              CALL MATVEC(WRK(IRES),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IW),1)
              CALL PRECONL(WRK(IW),WRK(IRES),IPAR)
          END IF

*  17. beta=||r||_2
          BETA = PDNRM2(LOCLEN,WRK(IRES),IPAR)

   20 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
      END IF

 9999 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE PIMDSETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,
     +                      NPROCS,PROCID,PRECONTYPE,STOPTYPE,MAXIT,
     +                      EPSILON)

c           PIM -- The Parallel Iterative Methods package
c           ---------------------------------------------
c
c                      Rudnei Dias da Cunha
c               Centro de Processamento de Dados,
c         Universidade Federal do Rio Grande do Sul, Brasil
c                              and
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c                          Tim Hopkins
c     Computing Laboratory, University of Kent at Canterbury, U.K.
c
c ----------------------------------------------------------------------
c
c  Description of parameter arrays
c   IPAR (INPUT)  : integer
c     ipar( 1): lda    (Leading dimension of a)
c           2 : n      (Number of rows/columns of a)
c           3 : blksz  (Size of block of data; used when data is
c                       partitioned using cyclic mode)
c           4 : loclen (Number of elements stored locally;
c                       *PARALLEL: Equal to at least m/nprocs or
c                                  n/procs depending if row or
c                                  column partitioning is used or,
c                                  in the case of cyclic partitioning,
c                                  it is a multiple of either
c                                  m/(blksz*nprocs) or n/(blksz*nprocs).
c                       *SEQUENTIAL: equal to n)
c           5 : basisdim (Dimension of orthogonal basis, used in
c                       GMRES)
c           6 : nprocs (Number of processors)
c           7 : procid (Processor identification)
c           8 : precontype (Type of preconditioning; one of
c                           0 : No preconditioning,
c                           1 : Left preconditioning,
c                           2 : Right preconditioning,
c                           3 : Symmetric preconditioning)
c           9 : stoptype (Type of stopping criteria used)
c          10 : maxit  (Maximum number of iterations allowed)
c
c   IPAR (OUTPUT) : integer
c     ipar(11): itno   (Number of iterations executed)
c          12 : status (On exit of iterative method, one of
c                        0: converged to solution
c                       -1: no convergence has been achieved
c                       -2: "soft"-breakdown, solution may have
c                           been found
c                       -3: "hard"-breakdown, no solution)
c                       -4: conflict in preconditioner and stopping
c                           criterion selected
c                       -5: error in stopping criterion 3, r^{T}z<0)
c          13 : steperr (If status is either -2 or -3, it gives
c                         the step number of the respective algorithm
c                         where a breakdown has occurred. Refer to the
c                         User's Guide for further information)
c
c   RPAR/DPAR (INPUT)  : real/double precision
c     dpar( 1): epsilon (The value of epsilon for use in the
c                        stopping criterion)
c
c   RPAR/DPAR (OUTPUT) : real/double precision
c     dpar( 2): exitnorm (The value of a norm of either the residual
c                         vector or the difference between two
c                         successive solution estimates according to
c                         the value of stoptype)
c           3,
c           4 : smallest and largest eigenvalues of Q1AQ2 (in the
c               symmetric case) OR smallest and largest real parts
c               (in the nonsymmetric case)
c           5,
c           6 : smallest and largest imaginary parts (only in the
c               nonsymmetric case)
c

C     .. Parameters ..
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=6)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION EPSILON
      INTEGER BASISDIM,BLKSZ,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE,
     +        PROCID,STOPTYPE
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION DPAR(DPARSIZ)
      INTEGER IPAR(IPARSIZ)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION EXITNORM
      INTEGER ITNO,STATUS,STEPERR
C     ..
      IPAR(1) = LDA
      IPAR(2) = N
      IPAR(3) = BLKSZ
      IPAR(4) = LOCLEN
      IPAR(5) = BASISDIM
      IPAR(6) = NPROCS
      IPAR(7) = PROCID
      IPAR(8) = PRECONTYPE
      IPAR(9) = STOPTYPE
      IPAR(10) = MAXIT
      IPAR(11) = -1
      IPAR(12) = -1
      IPAR(13) = -1

      DPAR(1) = EPSILON
      DPAR(2) = -ONE

      RETURN

      END
      SUBROUTINE PIMDTFQMR(X,B,WRK,IPAR,DPAR,MATVEC,PRECONL,PRECONR,
     +                     PDSUM,PDNRM2,PROGRESS)
      IMPLICIT NONE

*           PIM -- The Parallel Iterative Methods package
*           ---------------------------------------------
*
*                      Rudnei Dias da Cunha
*     National Supercomputing Centre and Mathematics Institute
*         Universidade Federal do Rio Grande do Sul, Brasil
*
*                          Tim Hopkins
*     Computing Laboratory, University of Kent at Canterbury, U.K.
*
* ----------------------------------------------------------------------
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
      INTEGER IPARSIZ
      PARAMETER (IPARSIZ=13)
      INTEGER DPARSIZ
      PARAMETER (DPARSIZ=2)
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),DPAR(DPARSIZ),WRK(*),X(*)
      INTEGER IPAR(IPARSIZ)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM2
      EXTERNAL PDNRM2
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONL,PRECONR,PROGRESS
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION ALPHA,BETA,C,EPSILON,ETA,ETA0,EXITNORM,KAPPA,RHO,
     +                 RHO0,RHSSTOP,SIGMA,TAU,TAU0,THETA,THETA0
      INTEGER BASISDIM,BLKSZ,CNVRTX,ID,IG,IH,IM,IM0,IP,IR,IRTILDE,ITNO,
     +        IV,IW,IXOLD,IY,IY0,IZ,LDA,LOCLEN,MAXIT,N,NPROCS,
     +        PRECONTYPE,PROCID,STATUS,STEPERR,STOPTYPE
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DOTS(1)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT,DSETRHSSTOP
      EXTERNAL DDOT,DSETRHSSTOP
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY,DINIT,PIMDGETPAR,STOPCRIT
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC DBLE,SQRT
*     ..

      CALL PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM,NPROCS,
     +                PROCID,PRECONTYPE,STOPTYPE,MAXIT,ITNO,STATUS,
     +                STEPERR,EPSILON,EXITNORM)

*  Check consistency of preconditioning and stop types
      IF (((PRECONTYPE.EQ.0).OR. (PRECONTYPE.EQ.2)) .AND.
     +    (STOPTYPE.EQ.6)) THEN
          ITNO = 0
          STATUS = -4
          STEPERR = 0
          GO TO 30

      END IF

*  Does not need conversion Y=Q2X for residual
      CNVRTX = 0

*  Set indices for mapping local vectors into wrk
      IR = 1
      IRTILDE = IR + LOCLEN
      IY = IRTILDE + LOCLEN
      IY0 = IY + LOCLEN
      IW = IY0 + LOCLEN
      IV = IW + LOCLEN
      ID = IV + LOCLEN
      IZ = ID + LOCLEN
      IP = IZ + LOCLEN
      IG = IP + LOCLEN
      IH = IG + LOCLEN
      IXOLD = IH + LOCLEN

*  Set rhs of stopping criteria
      RHSSTOP = DSETRHSSTOP(B,WRK(IR),EPSILON,IPAR,PRECONL,PDNRM2)

*  1. r=Q1(b-AQ2x)
      IF (STOPTYPE.NE.6) THEN
          IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IW),1,WRK(IZ),1)
              CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
              CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL DCOPY(LOCLEN,B,1,WRK(IP),1)
              CALL PRECONR(X,WRK(IW),IPAR)
              CALL MATVEC(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IP),1)
              CALL PRECONL(WRK(IP),WRK(IR),IPAR)
          END IF

      ELSE
*     r has been set to Qb in the call to dsetrhsstop
          IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
              CALL MATVEC(X,WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

          ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
              CALL PRECONR(X,WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IW),IPAR)
              CALL PRECONL(WRK(IW),WRK(IZ),IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)
          END IF

      END IF

*  2. w=y=r
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IW),1)
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IY),1)

*  3. v=g=Q1AQ2y
      IF (PRECONTYPE.EQ.0) THEN
          CALL MATVEC(WRK(IY),WRK(IG),IPAR)

      ELSE IF (PRECONTYPE.EQ.1) THEN
          CALL MATVEC(WRK(IY),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IG),IPAR)

      ELSE IF (PRECONTYPE.EQ.2) THEN
          CALL PRECONR(WRK(IY),WRK(IZ),IPAR)
          CALL MATVEC(WRK(IZ),WRK(IG),IPAR)

      ELSE IF (PRECONTYPE.EQ.3) THEN
          CALL PRECONR(WRK(IY),WRK(IG),IPAR)
          CALL MATVEC(WRK(IG),WRK(IZ),IPAR)
          CALL PRECONL(WRK(IZ),WRK(IG),IPAR)
      END IF

      CALL DCOPY(LOCLEN,WRK(IG),1,WRK(IV),1)

*  4. d=0
      CALL DINIT(LOCLEN,ZERO,WRK(ID),1)

*  5. tau=||r||2
      TAU = PDNRM2(LOCLEN,WRK(IR),IPAR)

*  6. theta=eta=0
      THETA = ZERO
      ETA = ZERO

*  7. rtilde=r
      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IRTILDE),1)

*  8. rho=dot(rtilde,r)
      DOTS(1) = DDOT(LOCLEN,WRK(IRTILDE),1,WRK(IR),1)
      CALL PDSUM(1,DOTS,IPAR)
      RHO = DOTS(1)

*  Loop
      STATUS = 0
      EXITNORM = -ONE
      STEPERR = -1
      IM0 = 1
      DO 20 ITNO = 1,MAXIT

*  9. sigma=dot(rtilde,v)
          DOTS(1) = DDOT(LOCLEN,WRK(IRTILDE),1,WRK(IV),1)
          CALL PDSUM(1,DOTS,IPAR)
          SIGMA = DOTS(1)

* 10. alpha=rho/sigma
          IF (SIGMA.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 10
              GO TO 30

          END IF

          ALPHA = RHO/SIGMA

* 11. y=y0-alpha*v
          CALL DCOPY(LOCLEN,WRK(IY),1,WRK(IY0),1)
          CALL DAXPY(LOCLEN,-ALPHA,WRK(IV),1,WRK(IY),1)

* 12. h=Q1AQ2y
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IY),WRK(IH),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(WRK(IY),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IH),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IY),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IH),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IY),WRK(IH),IPAR)
              CALL MATVEC(WRK(IH),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IH),IPAR)
          END IF

          DO 10 IM = IM0,IM0 + 1

* 13. w=w-alpha*Q1AQ2y0
              CALL DAXPY(LOCLEN,-ALPHA,WRK(IG),1,WRK(IW),1)

* 14. theta=||w||_2/tau0
              THETA0 = THETA
              TAU0 = TAU
              IF (TAU0.EQ.ZERO) THEN
                  STATUS = -3
                  STEPERR = 14
                  GO TO 30

              END IF

              THETA = PDNRM2(LOCLEN,WRK(IW),IPAR)/TAU

* 15. c=1/sqrt(1+theta^2)
              C = ONE/SQRT(ONE+THETA*THETA)

* 16. tau=tau0*theta*c
              TAU = TAU0*THETA*C

* 17. eta=(c^2)*alpha
              ETA0 = ETA
              ETA = C*C*ALPHA

* 18. d=y0+((theta0^2)*eta0/alpha)*d
              IF (ALPHA.EQ.ZERO) THEN
                  STATUS = -3
                  STEPERR = 18
                  GO TO 30

              END IF

              CALL DCOPY(LOCLEN,WRK(ID),1,WRK(IP),1)
              CALL DCOPY(LOCLEN,WRK(IY0),1,WRK(ID),1)
              CALL DAXPY(LOCLEN,THETA0*THETA0*ETA0/ALPHA,WRK(IP),1,
     +                   WRK(ID),1)

* 19. x=x+eta*d
              CALL DCOPY(LOCLEN,X,1,WRK(IXOLD),1)
              CALL DAXPY(LOCLEN,ETA,WRK(ID),1,X,1)

* 20. kappa=tau*sqrt(m+1)
              KAPPA = SQRT(DBLE(IM+1))*TAU

* 21. check stopping criterion
              IF (KAPPA.LT.EPSILON) THEN

*     r=Q1(b-AQ2x)
                  IF (PRECONTYPE.EQ.0) THEN
*     r=b-Ax
                      CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
                      CALL MATVEC(X,WRK(IP),IPAR)
                      CALL DAXPY(LOCLEN,-ONE,WRK(IP),1,WRK(IR),1)
                      CALL DCOPY(LOCLEN,WRK(IR),1,WRK(IZ),1)

                  ELSE IF (PRECONTYPE.EQ.1) THEN
*     r=Q1(b-Ax)
                      CALL DCOPY(LOCLEN,B,1,WRK(IZ),1)
                      CALL MATVEC(X,WRK(IP),IPAR)
                      CALL DAXPY(LOCLEN,-ONE,WRK(IP),1,WRK(IZ),1)
                      CALL PRECONL(WRK(IZ),WRK(IR),IPAR)

                  ELSE IF (PRECONTYPE.EQ.2) THEN
*     r=b-AQ2x
                      CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
                      CALL PRECONR(X,WRK(IP),IPAR)
                      CALL MATVEC(WRK(IP),WRK(IZ),IPAR)
                      CALL DAXPY(LOCLEN,-ONE,WRK(IZ),1,WRK(IR),1)

                  ELSE IF (PRECONTYPE.EQ.3) THEN
*     r=Q1(b-AQ2x)
                      CALL DCOPY(LOCLEN,B,1,WRK(IR),1)
                      CALL PRECONR(X,WRK(IZ),IPAR)
                      CALL MATVEC(WRK(IZ),WRK(IP),IPAR)
                      CALL DAXPY(LOCLEN,-ONE,WRK(IP),1,WRK(IR),1)
                      CALL PRECONL(WRK(IR),WRK(IZ),IPAR)
                      CALL DCOPY(LOCLEN,WRK(IZ),1,WRK(IR),1)
                  END IF

                  CALL STOPCRIT(B,WRK(IR),WRK(IZ),X,WRK(IXOLD),WRK(IP),
     +                          RHSSTOP,CNVRTX,EXITNORM,STATUS,IPAR,
     +                          MATVEC,MATVEC,PRECONR,PDSUM,PDNRM2)

*  Call monitoring routine
                  CALL PROGRESS(LOCLEN,ITNO,EXITNORM,X,WRK(IR),WRK(IZ))

                  IF (STATUS.EQ.-5) THEN
                      STEPERR = 21
                      GO TO 30
                  ELSE IF (STATUS.EQ.0) THEN
                      GO TO 30
                  END IF

              ELSE
*  Call monitoring routine
                  CALL PROGRESS(LOCLEN,ITNO,KAPPA,X,WRK(IR),WRK(IZ))

              END IF

*  y0=y
              CALL DCOPY(LOCLEN,WRK(IY),1,WRK(IY0),1)

*  g=h
              CALL DCOPY(LOCLEN,WRK(IH),1,WRK(IG),1)

   10     CONTINUE

* 22. rho=dot(rtilde,w)
          RHO0 = RHO
          DOTS(1) = DDOT(LOCLEN,WRK(IRTILDE),1,WRK(IW),1)
          CALL PDSUM(1,DOTS,IPAR)
          RHO = DOTS(1)

* 23. beta=rho/rho0
          IF (RHO0.EQ.ZERO) THEN
              STATUS = -3
              STEPERR = 23
              GO TO 30

          END IF

          BETA = RHO/RHO0

* 24. y=w+beta*y0
          CALL DCOPY(LOCLEN,WRK(IW),1,WRK(IY),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IY0),1,WRK(IY),1)

* 25. g=Q1AQ2y
          IF (PRECONTYPE.EQ.0) THEN
              CALL MATVEC(WRK(IY),WRK(IG),IPAR)

          ELSE IF (PRECONTYPE.EQ.1) THEN
              CALL MATVEC(WRK(IY),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IG),IPAR)

          ELSE IF (PRECONTYPE.EQ.2) THEN
              CALL PRECONR(WRK(IY),WRK(IZ),IPAR)
              CALL MATVEC(WRK(IZ),WRK(IG),IPAR)

          ELSE IF (PRECONTYPE.EQ.3) THEN
              CALL PRECONR(WRK(IY),WRK(IG),IPAR)
              CALL MATVEC(WRK(IG),WRK(IZ),IPAR)
              CALL PRECONL(WRK(IZ),WRK(IG),IPAR)
          END IF

* 26. v=Q1AQ2y+beta*(Q1AQ2y0+beta*v)
          CALL DCOPY(LOCLEN,WRK(IG),1,WRK(IZ),1)
          CALL DAXPY(LOCLEN,BETA,WRK(IH),1,WRK(IZ),1)
          CALL DAXPY(LOCLEN,BETA*BETA,WRK(IV),1,WRK(IZ),1)
          CALL DCOPY(LOCLEN,WRK(IZ),1,WRK(IV),1)

          IM0 = IM0 + 2

   20 CONTINUE

      IF (ITNO.GT.MAXIT) THEN
          STATUS = -1
          ITNO = MAXIT
          IF (EXITNORM.EQ. (-ONE)) THEN
              EXITNORM = KAPPA
          END IF

      END IF

   30 CONTINUE

      IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
          CALL DCOPY(LOCLEN,X,1,WRK(IZ),1)
          CALL PRECONR(WRK(IZ),X,IPAR)
      END IF

*  Set output parameters
      IPAR(11) = ITNO
      IPAR(12) = STATUS
      IPAR(13) = STEPERR
      DPAR(2) = EXITNORM

      RETURN

      END
      SUBROUTINE STOPCRIT(B,R,RTRUE,X,XOLD,WRK,RHSSTOP,CNVRTX,EXITNORM,
     +                    STATUS,IPAR,MATVEC,TMATVEC,PRECONR,PDSUM,
     +                    PDNRM)
      IMPLICIT NONE

*     .. Scalar Arguments ..
      DOUBLE PRECISION EXITNORM,RHSSTOP
      INTEGER CNVRTX,STATUS
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION B(*),R(*),RTRUE(*),WRK(*),X(*),XOLD(*)
      INTEGER IPAR(*)
*     ..
*     .. Function Arguments ..
      DOUBLE PRECISION PDNRM
      EXTERNAL PDNRM
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL MATVEC,PDSUM,PRECONR,TMATVEC
*     ..
*     .. Local Scalars ..
      INTEGER LOCLEN,PRECONTYPE,STOPTYPE
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION DOTS(1)
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT
      EXTERNAL DDOT
*     ..
*     .. External Subroutines ..
      EXTERNAL DAXPY,DCOPY
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC SQRT
*     ..
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DOUBLE PRECISION ONE
      PARAMETER (ONE=1.0D0)
*     ..
      LOCLEN = IPAR(4)
      PRECONTYPE = IPAR(8)
      STOPTYPE = IPAR(9)

      IF ((STOPTYPE.EQ.1) .OR. (STOPTYPE.EQ.2) .OR.
     +    (STOPTYPE.EQ.3)) THEN

*  Compute true residual if needed
          CALL DCOPY(LOCLEN,B,1,RTRUE,1)

          IF ((PRECONTYPE.EQ.2) .OR. (PRECONTYPE.EQ.3)) THEN
              CALL PRECONR(X,WRK,IPAR)
              IF (CNVRTX.EQ.1) THEN
                  CALL TMATVEC(WRK,XOLD,IPAR)
                  CALL MATVEC(XOLD,WRK,IPAR)
                  CALL DAXPY(LOCLEN,-ONE,WRK,1,RTRUE,1)
              ELSE
                  CALL MATVEC(WRK,XOLD,IPAR)
                  CALL DAXPY(LOCLEN,-ONE,XOLD,1,RTRUE,1)
              END IF
          ELSE IF (CNVRTX.EQ.1) THEN
              CALL TMATVEC(X,XOLD,IPAR)
              CALL MATVEC(XOLD,WRK,IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK,1,RTRUE,1)
          ELSE
              CALL MATVEC(X,WRK,IPAR)
              CALL DAXPY(LOCLEN,-ONE,WRK,1,RTRUE,1)
          END IF
      END IF

      IF ((STOPTYPE.EQ.1) .OR. (STOPTYPE.EQ.2)) THEN

*  ||r||<epsilon or ||r||<epsilon||b||
          EXITNORM = PDNRM(LOCLEN,RTRUE,IPAR)
          IF (EXITNORM.LT.RHSSTOP) THEN
              STATUS = 0
          ELSE
              STATUS = -99
          END IF

      ELSE IF (STOPTYPE.EQ.3) THEN

*  sqrt(r(Q1r))<epsilon||b||
          DOTS(1) = DDOT(LOCLEN,RTRUE,1,R,1)
          CALL PDSUM(1,DOTS(1),IPAR)
          IF (DOTS(1).LT.ZERO) THEN
              STATUS = -5
              RETURN

          END IF

          EXITNORM = SQRT(DOTS(1))
          IF (EXITNORM.LT.RHSSTOP) THEN
              STATUS = 0
          ELSE
              STATUS = -99
          END IF

      ELSE IF ((STOPTYPE.EQ.4) .OR. (STOPTYPE.EQ.5) .OR.
     +         (STOPTYPE.EQ.6)) THEN

*  ||Q1r||<epsilon or ||Q1r||<epsilon||b|| or ||Q1r||<epsilon||Q1b||
          EXITNORM = PDNRM(LOCLEN,R,IPAR)
          IF (EXITNORM.LT.RHSSTOP) THEN
              STATUS = 0
          ELSE
              STATUS = -99
          END IF

      ELSE IF (STOPTYPE.EQ.7) THEN

*  ||x-x0||<epsilon
          CALL DCOPY(LOCLEN,X,1,WRK,1)
          CALL DAXPY(LOCLEN,-ONE,XOLD,1,WRK,1)
          EXITNORM = PDNRM(LOCLEN,WRK,IPAR)
          IF (EXITNORM.LT.RHSSTOP) THEN
              STATUS = 0
          ELSE
              STATUS = -99
          END IF
      END IF

      RETURN

      END
